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ABSTRACT 


The minimum norm formalism of functional analysis is applied to 
the problem of controlling the neutron flux distribution in nuclear 


reactors. 


The application of the optimization technique is first described 
in terms of a general linear distributed parameter model. The general 
results thus obtained are later particularized to the specific case of 


a homogeneous slab reactor with the one-energy neutronic model, 


The reduction of distributed parameter systems into decoupled and 
more tractable subsystems is discussed. It is demonstrated how symmetry 
principles can facilitate the practical implementation of the optimization 
techniques to more realistic distributed models which may include for 


example reactor configurations in two or three dimensions. 


A new model expansion method is developed for solving the system 
of partial differential equations that describe the transient behavior of 
the most important fission product concentrations in thermal nuclear 


reactors. 
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CHAPTER I 


INTRODUCTION 


The thermal-power distribution in the core of a large nuclear 
reactor depends on the fission-reaction rate at each point of the core 


and it is neither homogeneous nor stationary. 


The rate at which fission reactions are produced is influenced by 
the neutron flux distribution. In turn this is affected by any change 
in the way neutrons are absorbed, diffused or produced throughout the 
core. These changes may be induced by a variety of causes. Power level 
adjustments, on-power refueling operations, fuel consumption and the 


accumulation of fission products in the nuclear fuel are typical examples. 


In most reactor designs, the neutron flux distribution is controlled 
by introducing (or extracting) a neutron-absorbent substance at specific 


locations of the core. 


The main objective for the reactor control system is to maintain 
the neutron-flux distribution within acceptable limits in order to avoid 


potentially harmful spots with high power density. 


It is recognized that optimal control theory can play an important 
role in improving the performance of existing reactor control systems 


and considerable research has been done in this area. 


The work presented in this thesis is a contribution to that research 


effort. 
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1 Scope of the thesis 


Chapter II is a review of the abstract minimum norm problems of 
functional analysis. The solution to particular cases which are utilized 


in later chapters are briefly discussed. 


Chapter III is devoted to the application of the minimum norm 
formalism to the problem of bringing the neutron flux distribution from 
an initial state into the vicinity of a desired distribution while 
minimizing a performance index that penalizes both the deviations of 
the flux from the desired distribution and the control effort. This 
problem is relevant to the particular situation where an external 
disturbance, such as on-power refueling, has affected the flux distribution 
and it is desired to reshape the resulting pattern into its original 
shape. This problem has been discussed in the past by several authors 
who utilized a variety of different optimization techniques. A review 
of their work is given in the chapter with the purpose of providing a 
yard-stick against which the optimization approach presented in this 


thesis can be compared. 


The application of the minimum norm approach is made first in 
terms of a very general linear distributed parameter model. The 
application is later particularized to the specific case of a homogeneous 
slab reactor with the one-energy neutronic model and a finite number 


of control rods. 


In Chapter IV several variations on the problem of controlling 
the neutron-flux are discussed. In addition to the requirements of the 
problem treated in Chapter III, the problem formulations in Chapter IV 


include new constraints imposed on both the total power generated in 


no! 


; are 
PA. ‘ on ot « 
Jo emsidorq sion muinbola sontteds S12 oe esas Bb vat co x9198 


ye 


=] 


bustii¢y or» fotviv Guado walooigng 63 wisn! th site! cong sn 


»fesr oath a aa ome sree 


Fae, ee 
asec mamtatm etia 2o worrasifeqe silty Basoweb et Tit joan 

et gotdgudiaviesl sul? movihen ads yatgrkad Jo peldosg oft oF 
Yo volatoly sH2 commit s383% “tw: int. 


~~ i> tia ke » 
Sizew GL 12 130222072 © nAVTIZSSD 5 i. 


= “a > 
¢ >. - 
> wee, 
| tao lsetveB edt dtodd esxliansqg Jets xebal sontari ots9G fy patie wheke 


" , e ‘ oi 1 - q an | oe 
ata? , mote Jortehs sfy bra ootsuddetats beviqab-sd9 aos} nef 


leeaeses ny ara poidansiea walisitseq sii of Jas sveleT all 


_ te 


jecdtyseth aul uta Heasostiga aod , gil Peaster “vewe Ti oral dows aoa 


toalcine 29 o¥ok weantssoa aptafmesz ef sascec7 z venta at 
440444 ° ‘= - » 1 »~ - 
—- - 

suoritus Iuteves Vd-gaeq 909 0) bosevsatb oved ead maidotg algt 


velver A ,esuntetets fnotderimtoqy toore?3th to “se! sav. & Wome 


seogzugq arid difw ‘tosqutin edd of aevig ef ow a 


e 
tp 
> 
my 
o 
— 
—_ 
&) 
be 
“~ 
“ 


2 tT a of, an 
PAS AL BEsneasig fanogcqas asolwesti2¢s3ns SAI A + rie Je riage 


a 


» DST RENO | baal 


ai 20322 afien 2) fotoseus aton sumiaim siJ3 Jo nok 38: ee 


otf .Jebue sedeumtaq besydixveit absall Leasaee’ vray #) 
1. ; lol 7 
Pty CD mg eT 5 Yo s@zo 3ittaeqe ssid 09 ae teal Ac heat ad; wile 


' + > "se an 
wo heen: ooluk? © fuse Toba oinesziier ¢gtus-6n9 sid EY 


the core and the state of the reactor. These constraints relate to 

the particular situation where it is desired to change or adjust the 

total power output from the core in a finite interval of time while 
minimizing both the distortion of the flux distribution and the control 
effort. It is demonstrated in this chapter how the additional constraints 


can be treated in the context of the minimum norm formalism. 


The application of the optimization technique is described in 
terms of a general linear distributed parameter model. Several examples 
are also given in which the particular case of a homogeneous slab reactor 


with the one-energy neutronic model is considered. 


Chapters V and VI are devoted to the development of computational 
methods that will facilitate the practical implementation of the 
optimization techniques to more realistic reactor models, which may 


involve for example geometrical configurations in two or three dimensions. 


In Chapter V the reduction of mathematical models into decoupled 
and more tractable subsystems is discussed. Thereduction procedure 


employed there relies on group theoretic arguments and symmetry principles. 


A short summary of the group theoretic results that are needed to 


describe the reduction approach is given for reference purposes. 


It is demonstrated in this chapter how by merely studying the 
geometrical symmetry of a nuclear reactor core it is possible to decouple 
the conditions for optimality that result from the application of the 


minimum norm technique. 


The reduction of a cylindrical reactor with the one-energy neutronic 


model is also given as an example. 
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A method is developed in Chapter VI for obtaining an approximate 
solution to a reactor core model that is suitable for describing the 


dynamic behavior of the reactor core during load-following operations. 


The method, which can be classified as a modal expansion approach, 
uses the eigenfunctions of the Laplacian operator and yields the 
solution to the model in an operator form that is amenable to application 


of the minimum norm technique. 


Both the symmetry reduction procedure of Chapter V and the modal 
expansion approach are combined in an example that involves a cylindrical 
reactor and utilizes the numerical data typical of a large, pressurised- 


water, natural-uranium nuclear reactor. 
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CHAPTER II 


THE ABSTRACT MINIMUM NORM PROBLEMS OF FUNCTIONAL ANALYSIS 


le Introduction 

In contrast with other optimization techniques, which generally 
approach problems of a particular dynamic nature, functional analysis 
with its geometric character provides a unified framework for optimi- 


zation problems of discrete, continuous, distributed or composite nature. 


This chapter constitutes a review of a class of functional 
analytical problems which form part of the well-developed optimization 
theory by vector space methods, the so-called abstract minimum norm 


problems. 


Although many optimization problems cannot be formulated in 
these terms, the availability of a large variety of different norms 
provides enough flexibility for the minimum norm formulation to be of 


importance, particularly in the field of optimal control. 


For reference purposes, the solution to certain particular cases 
which shall prove to be useful in later chapters, will be discussed in 


the next section. 
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29 The Class of Abstract Minimum Norm Problems 


Historically the method of moments may be seen to be the first 
functional analytical technique, used in connection with optimal control 
problems, in which the minimum norm formulation is employed. The origin 
of the method of moments can be traced back to the very early develop- 


ments in functional analysis. 


For reference purposes the problem of moments is stated below. 


Problem 1. The %-problem of moments 
Given the linearly independent elements fips ed cence EET 
a normed ltnear space B, the real numbers Oss 1=1,2,..+,N and the positive 
real number 2 find the necessary and suffietent conditions for there 


to extst a linear functtonal u defined on B satisfying the constraint 
u(f =O, Es hee fe ghey) (1) 


and wtth norm 


[ful] = eyp el <4 Q) 
tte] 
In his book Theorte des operattons Lineatres (reference [48], 
pp. 57, 74-75), Banach discussed this problem and mentioned the work 
of Riesz [63] and Helly [64] concerning the necessary and sufficient 
conditions for there to exist a functional satisfying the conditions of 


problem l. 


Krein [49] studied the problem in a very general manner and gave 


specific results for certain normed spaces. It is in his honour that 
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problem 1 is often referred to in this literature as "Krein's 


2£-problem of moments". 


The first application of Krein's results to optimum control theory 
were made by Krasovskii in connection with the time optimal control of 


a linear plant with amplitude-constrained control signals [51]. 


Krasovskii's ideas triggered the publication of a large number of 
papers which were concerned with amplitude, energy, power and fuel 
constraints on both the control signals and the system states as well. 
A modern exposition of the method of moments and/or a literature review 
on its application to optimum control theory can be found in the 


references [45, 53, 59]. 


For comparison purposes, the 2-problem of moments may be formulated 


in a different manner; 
Problem 1(a). 


Let x be a gtven vector of the n-dimensional eucltdean space 
EY let c be a convex set of linear functionals defined on the normed 
space B, 
|u(n) | 
e={u: |lu|| = su —— < 2, heB} (3) 
rn |[hl| 
then for a given linear bounded transformation F mapping ¢ into E, 
find the necessary and suffictent conditions for there to extst a 


uec satisfying 


Fu=x (4) 
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Motivated by the class of linear final value problems in 
optimal control, Balakrishnan discussed in his classical paper [52] an 


abstract minimum norm problem of a more general scope, 


Problem 2. 


Let F be a compact linear bounded transformation mapping a 


Hilbert space H, tnto another Hilbert space H Fora given. xc in Ho 


2° 2 


tt ts requtred to minimize 

[|Fu - ||? (5) 
subject to u betng tn the sphere ec in H, defined by 

[ul |? <2 (6) 


This formulation may be seen to be a variation on the 2-problem 
of moments in the sense that the control u is not required to be a 


solution to the equation 
Fu=x (7) 


but only the best approximation in the H,-norm. 


The compactness assumption on F gave this problem a very general 
scope since the range of F was no longer required to lie in a finite 
dimensional space. The solution to this problem was obtained by 


Balakrishnan [52] in the form given below. 
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Solutton to Problem 2. 
Etther 
sup | | te*rekr] 7 F*x|| < 2 (8) 
k>o 
tn whitch case the sequence 


-1 
Uy = [FPF W+kI] “Fx (9) 


ts such that Uy, converges to the optimal element u of minimal norm 


Lim ||Fu,-«||? = inf ||Fu-w||? = ||Fu || , (10) 
k>o uUee 2 

or 
sup | | (F*Peer] 77 Ftx|| > 2 (11) 
k>o 


in whteh case 
aly 
= x * 
u. [F F+k I] Fry ees: 


where k, ts adjusted so that [|u| |= ytelds the untque solution to 


problem 2, 


In a second paper Balakrishnan [58] extended these results to a 
Banach space setting and, using the theory of one-parameter semigroups, 
developed a general theory of optimum control for mathematical models 
described by the abstract Cauchy problem. Also in this paper, the 


compactness assumption on F was removed and, by invoking the properties 
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of convex functionals in Hilbert spaces, the solution to problem 2, 


given above, was shown to remain valid. 


Motivated by the basic mathematical structure of linear systems, 
Porter [60] presented a new approach to the general minimum energy problem 
in Hilbert spaces. Porter and Williams [55] extended the results of 
the previous paper to a Banach space setting and discussed the questions 
of existence, uniqueness and the properties of the optimal control. The 


problem treated by Porter and Williams is stated as follows; 
Problem 3. 


Let B and D be Banach spaces and T, a bounded ltnear transformatton 
defined on B wtth values in D, For each — tn the range of T find an 


element ueB that satisfies 


E=Tu (13) 


while mtntmizing the performance index 


| Joe | | (14) 


This formulation may be seen to be a generalization of the &-problem 
of moments in that the range of T is not required to lie in a finite 
dimensional space but in a general Banach space. The control u is 


required however, as in the problem of moments, to satisfy a given equation. 


In the earlier developments of Porter's work the assumption that 


T is a linear bounded and onto transformation is always present. 
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Knowing that distributed parameter systems have mathematical 
models which typically involve a transformation with dense range, 
Porter [57] removed the assumption of T being onto and considered the 
case of a linear bounded transformation with dense range. It is proved 
in his paper that the results previously obtained are still valid under 


the new assumption. 


Several variations on the general minimum energy problem, which 
dramatically enlarged the application scope of the abstract minimum normed 
formulations, were considered by Porter and Williams [56]. The most 


general of these variations is, for reference purposes given below; 


Problem 4. 


Let F be a bounded ltnear transformation from the Banach space 


B tnto the Banach space B,, let T be a bounded ltnear transformatton 


P 
from B onto (or with dense range) the Banach space D, and let t, ¥ and 


E be given vectors tn B, B, and D respectively. Find au in B satisfying 
Tu=t ele. 

that minimizes 
[[u-a| [74] [Fug |? (16) 


The questions of existence and uniqueness as well as the properties 


of the optimal control are discussed in references [56] and [32]. 


Of particular interest to the application considered in this thesis 


is the Hilbert space version of problem 4, in which case the unique 
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12 
solution obtained by Porter and Williams (see reference [32], pp. 348) 
is as follows; 


Solutton to Problem 4. 


The untque solutton u. of the Htlbert space verston of problem 


4 ts gtven by 


A x 
1 pt nek 9) (17) 


= 4P)~ 
u., [ I+F*F ] 
where n ts the untque vector in D sattsfying 


e=P[TeF*P] 2 (T  ns+P*y) (18) 


and where F* ts the adjotnt of F. T ts the pseudo inverse of T defined 


by 
Prear4[rT+] te (19) 


Provided that TT* ts tmverttble. 
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CHAPTER III 
OPTIMAL CONTROL OF DISTRIBUTED NUCLEAR REACTORS 
ee Introduction 


This chapter is devoted to the problem of controlling the neutron 
flux distribution in a nuclear reactor core where the spatial kinetic 


effects are important. 


The optimization technique of the minimum norm problem in Hilbert 
spaces is applied to the problem of adjusting the neutron flux for a 
general distributed nuclear reactor whose dynamic behavior is described 


in the neighborhood of an equilibrium condition. 


The problem consists of computing the control function that 
transfers the state of the system from an initial condition to a desired 
state, in a given period of time, and minimizes a quadratic performance 
index that penalizes the deviations from equilibrium so as to avoid 


potential spatial instabilities and high power density spots. 
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2 Background 


Kliger [65] in 1965 considered the problem of changing the neutron 
density level in such a way as to keep its spatial shape undistorted. 
The performance index to be minimized was the mean-square error in space 
and time, between the actual system response and the desired trajectory- 
The reactor core was modelled by the one-energy neutron group model with 
delayed neutron precursors. The model's parameters were assumed to be 
independent of the spatial coordinate,and a pseudocontrol function 
proportional to the product between the neutron multiplication factor and 
the neutron density was treated like a distributed control. Both the 
system state and the control were expanded in terms of an orthonormal 
set of spatial modes. With the spatial independence of the system 
parameters, a decoupled set of ordinary differential equations was obtained 
for the expansion coefficients. The control term for each one of the 
spatial modes was determined by making use of known results for point 
reactor models. The distributed pseudocontrol effect thus determined was 
later approximated by a finite number of spatially concentrated control 
rods. In spite of the simplicity of the model considered, Kliger's 
contribution holds the merit of having introduced several features which 
have been used by many authors, and in a sense became standard techniques, 
namely, the use of an orthonormal expansion and the treatment of a 


pseudocontrol function like a distributed input. 


In 1966 Wiberg [74,75,11] studied the optimal feedback control of 
spatial xenon oscillations. Wiberg considered a very general but linear 
reactor model and made use of a modal expansion in terms of the Kaplan 
modes [7]. The finality property that characterizes this set of modes 
allowed for the decoupling of the power control system from the spatial 


regulator. The spatially concentrated control functions were modelied by 
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means of a first order perturbation. Upon truncation to a finite 
number of modes and making use of a quadratic performance index, which 
penalized the state deviations from equilibrium and measured the control 
energy, the control problem was reduced to terms which are manageable 

by the theory of optimal control by state space methods [76]. Wiberg's 
Major contribution is the idea of decoupling the power control system 
from the spatial regulator-system. His chapter in [11] has become one 
of the most important references in the field. Hsu and Bailey [66] in 
1967 applied an extended version of the maximum principle to a simple one- 
group model for an homogeneous slab reactor with one delayed neutron 
precursor model. A quadratic cost functional was considered. Their 


approach is based on Wang's work on distributed parameter systems [46]. 


Stacey applied several optimization techniques to the problem of 
controlling xenon spatial transients in thermal reactor cores. In 1968 
[76], he treated the problem in terms of the dynamic programming 
formalism and gave a numerical example for a fairly realistic 3-dimensional 
reactor model. The large number of state variables present, however, 


made it necessary to severely limit the set of allowable controls. 


In 1969, Stacey [77] treated the xenon oscillation problem by 
using the calculus of variations for distributed parameter systems. In 
this case the optimality conditions were obtained in the form of a 
system of partial differential equations with mixed boundary conditions, 
which are very difficult to solve. A numerical example was given for 
the one-group model with xenon and iodine dynamics. A distributed control 
input and a quadratic cost functional were considered in the example, 


for which an approximate solution was obtained by applying the methods 
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of quasilinearization and nodal approximation. Knowing the large 
difference in the time constants of the neutron kinetics and the reactor 
poisoning process, Stacey assumed in his work that changes in the neutron 
flux and temperature distributions occur instantaneously. In 1970, 
Stacey [78] considered the application of variational synthesis to the 
optimal control of spatially dependent reactor models. The so-called 
method of variational synthesis consists of expanding the system state 
and the control input in terms of known functions of space and time. The 
conditions that would minimize a given performance index are obtained 

by simple calculus, in the form of algebraic equations. The choice of 
expansion functions in this case is arbitrary and the approach lacks the 


rigor of other methods. 


Chauduri [79] in 1972 applied the maximum principle to the problem 
of controlling xenon transients. A quadratic cost functional was 
considered and a distributed control input was assumed. Solutions were 
obtained for a one-dimensional reactor model via space and space-time 


discretization. 


Lazarevich [80] et. al. in 1972 considered the application of 
dynamic programming methods to the problem of controlling the flux 
distribution. Results were obtained for a slab reactor by applying 
Galerkin's method to the Hamilton-Jacobi canonical equations. An 
algorithm for computing the optimal control using the eigenvectors of 
the Hamilton-Jacobi operator was also proposed. The approach can be 


seen to be equivalent to Kaplan's modal expansion method yd ke 


In 1974, El-Bassioni and Poncelet [13] treated the minimum time 


control of spatial xenon oscillations via Pontryagin's minimum principle 
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and Kaplan's modal expansion method. The same problem was considered 
by Chiang et. al. in 1976 [15] but this time the lambda-mode analysis 


was applied in order to circumvent the problem of computing the Kaplan 


modes. 


In 1973 Iwazumi and Koga [18] considered the terminal control 
problem for a slab reactor with the one-group diffusion and one delayed 
neutron precursor model. The idea of introducing a pseudocontrol function 
was suggested, so as to replace the spatially dependent parameters of 
the model by both a combination of spatially independent parameters and 
a particular initial pseudocontrol distribution. The modal expansion 
method in terms of the "clean" reactor modes was applied and the familiar 
Riccati equations were derived for each mode. Motivated by the relatively 
large difference in the time constants of the prompt and delayed neutron 
kinetics, recently (1977), Assatani et. al. [19] suggested in connection 
with Iwazumi's approach, a singular perturbation method for solving the 


Riccati equations [82,83]. 


The first application of the abstract minimum norm formalism to the 
problem of controlling the neutron flux distribution was made in 1968 by 
Kyong [67]. Kyong treated the terminal control problem for a reactor 
core of cylindrical configuration, containing a finite number of control 
rods and with the neutron kinetics modelled by the one-group diffusion 
equation. Although the constraint on the maximum allowable control 
effort was not explicitly stated, Kyong's formulation of the problem may 
be seen to correspond to the abstract minimum norm problem treated by 
Balakrishnan in [52] and referred to as problem 2 in Chapter II of this 


thesis. The unique optimal control, for the terminal problem, was shown 
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to satisfy a coupled set of Fredholm's integral equations of the 

second kind. A rigorous method, based on the characteristic expansion 
[68], was proposed for solving the integral equations. It is clear 

from Kyong's work that the functional analytical formulation yields 
necessary and sufficient conditions for optimality in a form which is 
amenable to application of a different class of computational techniques, 
which in some cases may prove to be superior to the more conventional 


computational methods associated with the variational formulation. 


In 1973, Iwazumi and Koga also discussed the application of the 
minimum norm formulation to the terminal control problem, but this time 
a slab reactor model, with distributed control input and the one group - 


one delayed precursor neutron kinetics was considered. 


The relevance of the terminal control problem to reactor core 
control applications is questionable. Given that the neutron flux 
deviations from a desired state are only penalized at the end of the 
control interval, it is possible that large power transients could result 
from the implementation of the control function which is optimum for 
the terminal problem. This observation motivated the author to study 
the application of the minimum norm formulation to the case where the 
neutron flux distribution is brought from an initial state to a desired 
equilibrium distribution while penalizing both the deviations from 
equilibrium along the trajectory and the control effort. Some results 
of this investigation were published in 1977 by the author, Christensen 
and El-Hawary in reference [17]. This chapter and the next constitute 


an expanded version of that work. 
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as The Reactor Core Model 


Wiberg [11] described the mathematical model of a nuclear reactor 


core and discussed the nature of the most important kinetic processes. 


For convenience, the linearized model described by Wiberg will 
be considered here in its most general form. In order to avoid unnecessary 
repetition, the reader is referred to [11] for the details of the 
derivation. Also for reference purposes, the specialization of the model 
to the case of two neutron groups with delayed neutron precursors and 


xenon and iodine dynamics is considered in[7] and [11]. 


The dynamical equation of a general reactor near an equilibrium 


state is given by a linear vector partial differential equation of the 


form 


<— tr, PsA) eutr, cl e(rjauir) (1) 


where r is in general a 3-dimensional spatial variable, defined on the 
reactor extrapolated volume. w(r, t), for fixed r and t is an N- 
dimensional vector containing; the neutron flux at different energies, 
delayed neutron precursor densities, fission product and precursor 
densities, and fuel and moderator temperatures. A(r) is a matrix spatial 
operator involving the gradient and Laplacian operators. For fixed 

time t, U(t) isan M-dimensional vector representing the effect induced 

by the control rods. B(r) is a rectangular matrix spatial operator of 


the appropriate dimensions. 


Associated with equation (1) are boundary conditions of the form 
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pr, t)= 0 (2) 
at the reactor extrapolated boundary, and the initial condition 
W(r, ty)= Z (x) (3) 


at time ty: The neutron flux and current must also satisfy continuity 


conditions at internal boundaries. 


In general,one wishes to find a solution to equation (1) ina 
normed function space of interest. Here the particular case is considered 
where for each fixed time t, the state W(r, t) is a function in a real 


Hilbert space H, endowed with an inner product of the form 
f a. 
< Z(t), W(t) >, = J Wr, t) Q(r) 2, t)dr (4) 
V 


in which the integral is over the core extrapolated volume V. Q(r) 
is an N-dimensional, positive definite matrix with space-dependent 


: ik 
entries. W denotes the transpose of W. 


Similarly the control function U(t), for a fixed time t, is 


assumed to be an element of a finite dimensional, real Hilbert space E, 


with inner product 
At 
EPG oe® ae, oo Y Ce) R UC) (5) 


where ye is the transpose of Y, and R is a positive definite matrix of 
the appropriate dimension. The control functions considered belong to 


the class of functions with a bounded norm of the form 
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f wlce) R v(e)dt = |Juce)| | ,at (6) 


ct 
oO = 
fh 


The assumption is made that equation (1), together with the 
associated boundary conditions, constitute a well posed abstract Cauchy 
problem in relation to the norm topologies of H and E. Therefore the 


unique solution to (1) and (2) can be written in the form 


c 
Berea y=c(rsitt r!, th) Z(ni} + | G(r, ts c', 1) BG) Ulo)dr 


t 
: (7) 


where G(r, t; r', t) is the strongly continuous semigroup associated with 


the operator A(r), which satisfies, 
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eitemts rit) a=. Identity (9) 


and is also a strong solution to 


<— G(r, Pir itac yw (ria teeter fy Acuna) 


(10) 


in the Hilbert space H. 


Equation (7) is the solution to (1) in the sense that yp satisfies 


the relations 
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, 1 
lim ||h “[y(t+h) - v(t) ]- A v(t) - B u(t) ||, = 0 (11) 
h>0 
+ 
formed it > to» and where h > O01 denotes h > 0 through positive values, 


and 


lim ||y(t) - z.||. = 0 (12) 
£0, o''H 


for all initial conditions Zo in the domain of A; 


D[A(r)] = {Z(r) 3 | Jar) 2(r)| |, Sim gandazinpa= 0 
in the boundary of v}. 


Several remarks are pertinent; 
1. - Balakrishnan [58] and Wang [41] discuss the well-posedness 
question in connection with the abstract Cauchy problem in Banach spaces, 
and give sufficient conditions for the unique solution to be in the 
form of equation (7). Namely it is required that: a) A(r) is the 
infinitesimal generator of the semigroup G, b) BU is strongly measurable 
and Bochner integrable in every finite time interval, c) BU is in 
D[A] for almost every t and d) ||ABU|| is integrable in each finite time 


interval. 


A detailed discussion on these definitions and properties can be 
found in reference [89]. 
2. - When the matrices Q and R, appearing in the inner product 
definitions (5) and (4), are chosen to be the identity in their correspond- 


ing spaces, it is clear that the solution (7) would satisfy (1) in the 
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mean square sense and the control functions being considered would belong 
to the class of controls with limited energy. The inclusion of a matrix 

Q other than the identity, provides a wider range of norms, which would 
allow for example, to emphasize the difference between the state 
distributions and the equilibrium condition in a given region of the 
reactor where core-damage due to large flux deviations is more likely 

to occint 

3. - Although it is convenient to describe the application of optimization 
techniques to reactor control problems in terms of a general model, that 
at least in theory applies to all possible situations, it should be pointed 
out that a model of such generality is seldom used in reactor simulations, 
for the simple reason that it is too complex to treat computationaly. 

In practical situations, approximations of diverse nature are always 
introduced. Of these, the most common is the hierarchical decomposition 
of the general model according to the time scale that is relevant to the 
particular situation or problem in case. 

4. - If the conditions of remark 1 are satisfied, then the solution to 

the mathematical model exists and it can always be obtained in the form 
of an integral expression of convolution type such as (7). However, it 
should be pointed out that the conditions of remark 1 are difficult to 
verify in practical situations. Also, obtaining the solution in the form 
of (7) is seldom an easy task. 

5. - Finally a few explanatory remarks on notation are mandatory. The 
semigroup G(r, t; r! t) represents in our case an integral operator. The 


t 


kernel of this operator depends on the four variables r, t, r' and tT. 


Of these, only r' defines the dummy variable of integration. 
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4. The Optimal Control Problem 


This section is devoted to the control problem where the state 
distribution is brought from an initial state to a desired equilibrium 
condition, while penalizing both the deviations from equilibrium along 


the trajectory and the required control effort. 
The performance index considered here is given by 


t 


1 
R= ellyirs ey — Zr) IP OG luis tye Cr) lardt 
t V 
0 
t 
1 it 
+ f  U-(t) R U(t)dt (14) 


where Z is the desired state. Ko is a positive constant chosen so as 
to establish the relative weight of the second term in (14). w is the 


system's state given by equation (7). Q and R are defined as before. 


It is convenient to introduce the following notation: H, and H 


uf 2 
will denote the Hilbert spaces with inner products 
at 
Ri iy>ree pyer<U(c pay(t)> dt (15) 
ik ihe 
0 
and 
ih 
<Z, Wo, = if rage W(t)>,,dt (16) 
2 to 


respectively. Also, the operator F, mapping Hy into Hy will be defined 


by 


ro 
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t 
F(r, t; t) U(t)= f G(r, t3 x', tT) Br") U(t)dt (17) 
t 
0 


With this notation, the cost functional given in (14) can be re- 


written in the form 
J(u) = ||FU-Z+ CyZo| hn, + ky | IU rs (18) 


where Zo is the state distribution at time to and Gy is the semigroup 


defined by 
= ° ' 
Gy Ciresteare ty) (19) 


which maps the space Dato itself. 


The minimum norm formulation of the control problem can be seen to 
correspond with problem 2 of the previous chapter. Therefore the optimal 


control U can be determined from 


[F*F + ko I]U = F*(Z - C529] (20) 


where F* is the operator adjoint of F. 
4,1 The Adjoint F* 


The adjoint F* is related to F through the inner product relation 


Bee <U, iE, (21) 


the left side of (21) can be expanded as follows 
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1 t 
f < f Ces ens (t) dt. Z2(t)>,dt (22) 


Cen at) 


where for simplicity only the temporal dependency is explicitly shown. 
The inner product in H involves an integral operation throughout the 
reactor core volume. After interchanging the order of integration, 


one finds that (22).can. be written in the form 


E, 

Lac 
f f <Gits tT) BUCT), Z(t)>,drdt (23) 
X9 “9 


which in terms of the adjoint operators G* and B* becomes 
ct a 


SEU, 22 c= f <U(t), BG*(t3 Tt) Z(t)>,dtdt (24) 
t 


0 to 


and after interchaning the order of integration 


hg Tass 
SFU 82>. ff <U(t), B&G*(t; tT) Z(t)>,dtdt (25) 
2 t. T 
0 
ay ea) 
= f <U(t), f  BRGR(tst) Z(t)dt>,dt (26) 
t . 
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from which F* is obtained as follows, 


sh) 


F*(t, t)Z(t) = f B*G*(t; t) Z(t)dt 
44 


Similarly, B* is related to B through the relation 


<BU, Z>_= <U, BRZ> 


H 


Expanding the left sideof (28) one finds that, 


<BU, Z>.= f 2'(r) Q(x) B(r)Udv 
V 


§ tplcey air) z(r) J av u 
V 


rR? f Ble) Qhir) z2(r)dv) "RU 
V 


= ak 
<u, Rf B(r) Q’(r) Z(r)dv>, 
V 
: * . . 
from which it follows that the adjoint B is given by 
Bits, Rar ef Bi Cr) oh (r) Z(r)dv 
V 


where R has been assumed symmetric. 


(27) 


(28) 


(29) 


(30) 


(31) 


(32) 


(33) 
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4.2 The Necessary and Sufficient Condition for Optimality: 


A Fredholm's Integral Equation 
Substituting F and F*, as obtained in expressions(17) and (27), 
x 
one finds that the term F FU can be expanded in the form 
rl 


t 
FAFU = f B*G(t; t) f G(t; a) BU(a)dadt (34) 


T Co 


where again for simplicity only the temporal dependency is explicitly 


shown. 
Interchanging the order of integration in the expression (34), one obtains 


c 


T 1 
F*FU = f f B*G*(t; tT) G(t; a)Bdt U(a)da 
tout 
0 
Seert 
Ss f BAG*(ts tT) G(t3 a) Bdt U(a)da (35) 
gy a 


which, in turn, can be rewritten in the form 


t 
ik 
F*FU = f K(t; a) U(a)da (36) 
0 


where 


t 
f 1 BYG*(t, peo eee bat tor d<T 
T 
K(t, a) = (37) 


t 
f 3 BXG*(t; Tt) G(t; a)Bdt for art 
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It follows from expressions (20) and (36) that the condition for 


optimality is given by a Fredholm's integral equation of the second 


kind, 


t 
1 
U(r) = A(t) - =f K (ts a) U(a)da (38) 


0 to 


where the function A is defined by 


uae 
T™) = k [Z - GoZ9] 
0 
aT 
Ge ech Cc: a 
Ey ey BG (t3 t) [Z- G)Z)]dt (39) 


It can be noticed from the defining relation (37) that K(t, a) is 
not in the form of a product or a sum of products of functions that 
depend on a single variable. In other words, K(t, a) is a non-degenerate 


kernel. 


It is worth mentioning here that the necessary and sufficient conditions 
for optimality derived by Kyong [67] for the terminal control problem 


are in the form of an integral equation with a degenerate kernel. 


Since the key ingredient in the characteristic expansion method 


proposed by Kyong for solving the optimality conditions is precisely 
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the degeneracy of the kernel, it follows that this approach is not 


applicable to equation (38). 


It should be pointed out however, that equation (38) is amenable 
to application of approximating techniques such as contraction mapping 
algorithms [90], or methods that reduce the integral equation into 
a system of algebraic equations via temporal discretization or function 


expansion techniques [27]. 


4.3 An Example 


Consider an homogeneous slab reactor model and the one-neutron 
group diffusion equation 
fs M 
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with boundary condition 
PLOs7C) <= Ob, £)*="0 
and the initial neutron flux distribution $¢(r, to) > 


D is the neutron diffusion coefficient, Le and te denote the fission 
and absorption macroscopic cross sections respectively. V is the average 
neutron velocity, v denotes the number of neutrons generated per nuclear 
fission. b denotes the reactor width. The control absorption cross 
section at the location tr; is represented by u, (t). 
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penalizes the flux deviations from the equilibrium distribution o, (r) 
6 (r)= oly = sin@] (42) 
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where ou is a constant. Equation (42) its a solution to 
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0 =D a ea [v2_- Z1¢ () (43) 


subject to the boundary condition. 


Small flux deviations from o, (x) 


Virg t)=.¢6(r, -t)- oe) (44) 


satisfy the linearized equation 
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which can be solved by applying the method of separation of variables. 


The solution is given by 


b ' (Lat) 
v(r, t)= L = sin( “p OLS sin(— ) Z)(r')dr'Je ms S 
£ A (t-t) M nmr, uo 
+ (21°! v6, SL Let “sin@Ft) since sing Su, (16s 
0 n i=l 
(46) 
where 
Aas (vip- 2,- ED? div (47) 


In terms of the notation introduced at the beginning of this 
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ine operator G(r, ts xr", ty)» in this case, is self adjoint. 


The adjoint F* is given by anM-dimensional vector operator with entries 


£*(T5 r, t) defined by 
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d (2, -tT-a) d toa | 
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where |-| denotes absolute value. 


The function A(t) defined in (39) becomes in this case, 
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and making use of (56) and (49), the ith entry A, (t) is obtained as 


follows 
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and becomes, after integration, 


A, (t)= 
, i age s 68 
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The optimality conditions are finally obtained in the form 
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4.3.1 Discussion 

More conventional techniques [11] approach this control problem 
through modal expansion methods and invoke variational principles which 
yield necessary conditions for optimality in the form of an infinite 
system of ordinary differential equations with mixed boundary conditions. 
The question of how many modes should be included in a finite dimensional 
version of the system is well known to constitute a very difficult problen, 


which in most cases can only be solved through trial and error. In contrast, 
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condition (63) forms a system of M integral equations with non-degenerate 
kernels. Although these kernels are represented by infinite series, 

the question of how many terms should be considered in a finite series is 
a tractable problem and, as it is shown below, error bounds can be 


estimated with relative ease. 


Consider the kernel 
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which is obtained by truncating the series after the first P terms. 


The error is given by 
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is finally obtained. This estimate is a function of P, which rep- 


resents the number of terms considered in the finite series approximation 


to Sie a). 


4.3.2 A Method for Computing the Optimal Control 


Equation (63) can be rewritten in a more compact form 
U(t) = L(t; a) U(a) (76) 


where U(t) is an element of the space H, and L is a nonlinear trans- 


1 
formation from Hy into itself defined by 
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L(t} a) Ula) =A(t) - = f k(t; a) U(a) da (77) 
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and A(t) is a function in Hy with entries A, (7) defined by equation (62). 
K(t; a) is a matrix operator with entries Sine a) defined by equation 


(64). 


Invoking the contraction mapping theorem [90], it follows that 


the sequence 
u(t) = L(t; a) U 5-1 (78) 


based on the initial guess Uy (t) in Hy would converge to the unique 
solution to equation (76) provided that L is a contraction, that is, for 


any two functions U and W in H> 
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where C is a constant and 0<C<l. 


In what follows, a rough estimate of the constant C will be 
computed in order to test for the contraction condition and establish 
a criterion for the applicability of the successive approximation 


algorithm (78) to equation (76). 


Consider the norm 


t 


: Z 
| | f K(t3 a) U(a)da| |, (80) 
1 


=n 
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It is clear from (93) that a sufficiently large value of 


Ko would make C less than one. 


The design parameter ky? which was introduced in the performance 
index (41), provides the flexibility required to obtain the desired 
balance between system response and control action. A suitable value 
of ky) can only be found through systematic searching which of 


necessity involves repeated computations for different values of be 


Making the supposition that the contribution of the flux 
deviation term at its maximum in the performance index should equal 
the contribution of the control energy term at its maximum, Wiberg [11] 
has outlined a method for obtaining an initial rough estimate of ke 
For the case where the flux deviation is to be no more than 10% of the 
steady state flux and where the control reactivity is to be no more than 
10% of the total material reactivity, Wiberg proposed the initial rough 
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where M is the number of controls and the other variables are defined 


as before. 
Substituting (94) in (93) the following result is obtained, 


4 (t,-t.)b 


fe =ta} 120 
10° 4 1.97 x 1072 


2 
C = = [(vE,-2 VIE 
(95) 


Although (95) is a very conservative estimate, it provides a useful 
criterion to determine whether or not the contraction mapping algorithm 
would yield a convergent sequence. If C is less than one and U denotes 
the solution to equation (76) then the contraction mapping theorem 


ensures that the sequence 
U =LuwU (96) 
would converge, with the rate of convergence given by 


| |u-u_ || [LU ,-U ol ly (97) 


n 
a a 1 


where U is the initial guess. Consider for example the reactor 
Oo 
-1 
constants given by Wiberg (page 359, reference ri}; (vz .-2)V=0.256sec x 
b=250cm, M=2 and VD=1600cm”/sec. In this case, the criterion (95) 


would take the form 
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aa - 2 
¢ = 5.26x107*[0.5 [t,-ty]* + 6-11x107 [t,-ty]]” (98) 
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and since the time scale of the nuclear model considered is of the 
order of seconds, it can be concluded that the successive approximation 


algorithm (96) is applicable to problems in this range. 
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CHAPTER IV 


FURTHER APPLICATIONS OF THE MINIMUM NORM FORMULATION TO PROBLEMS IN CONTROL 


OF DISTRIBUTED REACTORS. 
ims Introduction 

Several interesting variations on the optimal control problem of 
Chapter III would result if in addition to the requirement that the control 
function minimizes a performance index, which takes into account the 
deviations of the system's state trajectories from a desired distribution, 
the state of the reactor core is also required to satisfy one of the 
following conditions: 

(a) At the end of a given time-interval the total power 


output from the reactor core should be equal to a specified value. 


This constraint is relevant to the problem of adjusting the power 
level while minimizing both the control effort and the distortion of the 
flux distribution. 

(b) The total power output from the reactor core should match 


a specified load trajectory along a given time-interval. 


This constraint is relevant to the problem of controlling a nuclear 
reactor core during load-following operations. 
(c) At the end of the specified time-interval the state of the 


system should match a given distribution. 


Although this constraint is also relevant to the problem of adjusting 
the power output from the core, it is more restrictive than case (a), 
stated above. 


The present chapter is devoted to the application of the minimum 


norm formulation to these particular problems. The scope of the chapter 
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is focussed upon the derivation of necessary and sufficient conditions 


for optimality in a form that lends itself to computation. 
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tds The Nuclear Reactor Model 


Although an accurate description of the power generated in a nuclear 
reactor core would have to take into account both the continuous energy 
spectra of the neutrons in the core and the energy-dependent fission cross 
section of the fuel, adequate results are often obtained through the 
simplifying assumption that the neutron population can be grouped into a 
number of monoenergetic groups whose dynamics are described by multigroup 


diffusion theory [37]. 


Here it is assumed that the total power produced in a nuclear reactor 


core is given by 


p(t) = 


hos Zz 
oq 


eermen es (ele (r,t) de (1) 
Voki 4 


where the integral is over the volume of the reactor core and ep is the 
i 
energy released per fission involving a neutron in the ith-energy flux 


.- ue is the fission macroscopic cross section for neutrons in the ith- 
a 
energy interval. ne denotes the number of energy intervals which cover 


the neutron energy range. 


This assumption is often employed in engineering applications. The 


reader is referred to Chapter IV of reference [91] for a detailed discussion 


of the one-energy group approximation. 


As was the case in Chapter III of this thesis, here the assumption 
is made that near an equilibrium condition the state, w(r, t),of the 


nuclear reactor core is given by 
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where all the variables are defined as before. 


In terms of the state wp(r, t) the total power p(t) can be modelled 


by the expression 


p(t) = f H(r) vr, t)ar (3) 
V 


where H(r) is a space-dependent row matrix with the appropriate dimension. 


In what follows it is assumed that the state wW(r, t) and the control 


U(t) belong to the Hilbert spaces H and Hy respectively. Hy and H, are 


endowed with the inner products 


t 
1 
Velen) <V CE), Z(t)> de (4) 
hss oats 
and 
aT 
aly Weare f k<U(t), W(t)>,dt (5) 
1 to 


where ko is a positive real number. The inner products in H and E are 


defined by 


<v(t), Z(t)>_ = f 2 (r, t) Qe) vr, t)dr (6) 
V 
and 
<U(t), W(t)>, = W(t) R U(t) (7) 
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in which R is a M-dimensional positive definite matrix. 


space-dependent, N-dimensional positive definite matrix. 
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at Optimal Control of the State Distribution with Power Level Adjustment 


Suppose that at time to the total power output corresponding to the 
initial state Z (x) is Po: The problem of finding the control function 
U(t) that changes the total power output from Po to Py in a given time 


interval [tys ti], and minimizes the performance index 


t 


1 
su) =f f (wr, t) - zr, t)]? ae) [v(r, t) - Z(r, t)]drdt 
to V 
t 
l T 
+k, f U (t) R U(t)dt (8) 
t 
0 


where Z(r, t) is the desired state distribution, will be considered in this 


section. 


It should be noticed that although no explicit constraint is imposed 


on the total power output 


p(t) = J H(r) Clr Cees ty) Z)(r')dr 


y. 
Pole ede) Greets rue eat UCT ydrdr (9) 
Vaeet 


along the given time interval, implicitly, the performance index (8) does 


penalize the output power variations along the trajectory. 


The present problem includes on an equal footing the problem of 
controlling flux oscillations, minor power adjustments and load following 
operations. The particularization of the problem into any of these cases 


depends on both the time scale and the reactor model considered. 
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In order to be able to profit from the results discussed in Chapter 
II, one should consider first the formulation of the problem in the 
context of functional anaylsis. With this in mind the following definitions 
are introduced: F denotes the transformation from H, into H., defined by 


u 


t 
F(r, t3;t) U(t) = f G(r git Sear eo re (Td Tee (10) 


C9 


T denotes the transformation 


Le A 
Pttwa tuo), = f f H(r)-CUrar er mptrs), UCo}drdr, 
b V 
0 
(11) 
Also, y(r, t) and €(t) are functions defined by 
Uirentoeses (Ceetoe—aC(rsits ri, ty) Zr") (12) 
and 
E(t) = p(t) - [ H(r) G(r, ts r', ty) Zy(r")ar (13) 
V 


It should be noticed that at time ty the function E(t,) is completely 
specified. The transformation T and the function & when evaluated at time 
t. will be denoted by qT, and ei Finally, H., represents the one-dimensional 
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real Hilbert space endowed with the inner product 
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With these definitions the optimal control problem of this section 


can be formulated as follows: 


Find the control U in Hy that minimizes 


50) = |lru-yll, + [lvll, (as) 
2 1 
and satisfies 
oi = TU (16) 


where T) is a linear and bounded transformation from Hy onto H,. 


This formulation is recognized at once as Porter's abstract 


minimum norm problem. 


SHANE The Necessary and Sufficient Conditions of Optimality 


By invoking Porter's result it follows that the control U which 


minimizes the performance index (15) and satisfies (16) is given by 
-1 ,+ 
U = [1+ F*F] (Tjn + F*y) (17) 
where n is the unique element of H., satisfying 


a Uk (18) 


+ 


The adjoint F* transforms H, into Hy: The pseudoinverse Th> 


associated with Ty» is given by 
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+ x =] 
T, = T,[T,T,] (19) 
} k 
and the adjoint qT is a transformation which maps A, into Hy: 


Since there is no straightforward approach to computing the inverse of 
x 
[I + F F] it is necessary to express conditions (17) and (18) in a more 
standard form, amenable to computation. To this end, condition (17) is 


rewritten as follows 
k ae k 
U=-F FU Tin ce ty aye (20) 
Also, given the algebraic property of the pseudoinverse 
+ 
iT WN aka (21) 
and making use of (18), one finds from (20) that 
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1 Hi 


Solving for n from (22) and substituting it in (20), one obtains 
the following operator equation 
= trent rou +3 A (23) 


+ 
U = [T, T, 


in which the function A is an element of Hy defined by 
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- s = 24 
f= 18 > (2, 7, 021 Fy (24) 
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The reader may notice that (23) lends itself to application of 


successive approximation techniques. It is true that the pseudoinverse 


a : : : ° ? % 

qT, implicitly involves the inverse of the operator TT: However, given 
* 

the fact that qT, has finite dimensional range, it follows that T,T, is an 


operator with finite dimensional domain and range. More explicitly, 
* 
TT, transforms the one-dimensional Hilbert space H, into itself. The 


* =] 
inverse [T, T)] » therefore, should pose no computational difficulties. 


In order to fully characterize the optimality condition (23), it is 
k 
still necessary to derive explicit expressions for the adjoint F and the 


+ 
pseudoinverse T , 


* 
a2 The Adjoint F 


Following the procedure outlined in section 4.1 of Chapter III, 


* 
it can be shown that F is given by 


1 
* 2h 
Hts (ton t EZ (tC) mak, CEB Gite eZ (ede (25) 
T 
where for simplicity only the temporal variable is explicitly shown. The 


* 2 
adjoint B transforms H into E according to 


Bz =R ! f Bl(r) Q’(r) Z(r)dr (26) 
4% 


* : 
and the adjoint G transforms H into itself. 


Replacing Z(t) in equation (25) by the defining relation (10), one 
obtains 
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1 
FY (cs t) F(t3 a) Ula) =k f Be (t; t) f G(t3 a) B U(a)dadt 
4f | 
0 
(27) 


to gobdesdiigae 02 hraadt abunal (OS) ome ; 


Sprsvilabysey, S07 Jens sites atl) .weuphedggl aolsnety wag 
| ES | 

6 cat a a 

evig <tVSOH a TT sczareqo sits Se covet ae vovio 
oo wt TT sede @Woblot a! .vetms be: ro nani dak end pe 4 
: «& ~ a, y 
btobotieve w¥hk° .epeny Ere noaedion far athe sitet « 104 


47 - .3fenet o9G! .H sowre Syacdlih ienalseaoes me v9 sis varta om 


: at 5 yy 
a@ekt rorvikh Lesolsaivuiroivo on veod blwars atoms rac . 


i hi 2, 
at 32 . (22) qnidthnes 22°) om) 190 oft aalroteetado Vile? of 19 


pid fies Intetha os ore ..que ovigst a4 ys 
as 
- Taf 
- 2 a a. 
1 gesotbA off See 
* > 
os | of to 2.8 totons i os leg ry ba org | odd pulbwod 


~~ 
ae * 


. 


oeviy at TI tac? ougied 


oa 
o*4. 
eo. 

au ee 


oO 
) 


(28) ib(ais (r wa" va (s)8.(a oe 
ee 


7 , , ‘ 
- 7 ay 
+ ue 
°@ . iM ‘a fy _ 

fi ,uwode yiottottoxs pr? eicecsav Lerog@e) aot yigo ¢2 to tiga S. 
es if 

x iJ as ‘ 

~ a 

. 


‘$ getbvesse 2 ovat H aernpoctesai aa 


(a5) rb{3) tay" ¢ 


iteestt opal a aan “02 


<3 


“a ab bas ? (#8) noes a 


- BRO AOR LQOnrnie Kee 
« - ed i 


4 


where again for simplicity only the temporal dependency is explicitly shown. 


After interchanging the order of integration in expression (27), 


it follows that 


t 
* ha 1 
Fo(t3 t) F(t; a) U(a)=ko” f | K(t3 a) U(a)da (28) 
(e 
0 
where 
t 
1 x x 
f BG (t; t) G(t3 a) Bdt for a<t (29) 
K(t3 a) = 
t 
1, x 
f BOG ats 97) 0G (tse ojepdL Lor dot (30) 
fo} 
2 of 
Sys. The Pseudoinverse Tj 
x 
The adjoint Ty is related to T, through the inner product relation 
* 
<T, Us age = <1) T) ane ‘ (31) 


Using the expression for T, given by (11), the left side of (31) 


1 


is expanded as follows 


f f H(r) G(r, ty3 Te CSO rae Corrie s (32) 
0 V 


which, in terms of the inner product in H, takes the form 


2 -1°T 
i <G(t,3 tT) BU(t), Q HE>, dt. (33) 
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From (33) it follows that 
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M -1,* * -1°T 
<T, y, sik = <U, Ko BGQ éH otis . (34) 
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The adjoint qT) 


is therefore given by 


k -1 ** = 
T, (t) =k : BG (t)3 TO “ 


oath 
0 H 


(35) 


x 
Upon substitution of the expression for B , given in (26), (35) leads to 


* nTsP— k a a 
Teast Ry (0B (x) Q(z) Gi@, ty:kles) Ox (cH lade 
V 
(36) 
* 
In view of (36), the term TT, becomes 
* ~1 
oa eo (37) 
where d is the real number defined by 
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ie V V 
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be 1 =e ats aa t! 
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* Rone 
It follows from (36) and (37) that the pseudoinverse T,{T,T,] 


can be finally obtained in the form 


T(t) = te K, ) (39) 
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where K (t) is an M-dimensional column matrix defined by 


K(r) = 8 f BG) Qh) CG, ty x1, 1) Ce") Hear 
V 


(40) 


ag The Optimal Control: The Solution to a Fredholm's 
Integral Equation 
* + 
With the help of the explicit relations for F and Ty» it is now 


possible to fully characterize the optimality condition (23). 


Recalling expression (11) for T, and expression (28) for F*FU 


one obtains 


a 


ru = koe 
ib See 3 K,7 (a) Ula) d(a) (41) 
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where the kernel Ry ™ is defined by 
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t V 
0 
(42) 
and, in view of (39), it follows that 
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Ys eee : Retr eotay Ufo). de (43) 
T, (t) TF FU = (d Ky) ‘ - T II : 
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Finally, making use of (43) and (28), one finds that the optimality 
condition can be reformulated in terms of a nonhomogeneous Fredholm's 


integral equation of the second kind, 
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U(t) = A(t) + ko f [d °K, (r) K, 7 (a) - K(t; a)] U(a) da 


(44) 


where A(t) is the function defined by (24). 


It would be interesting to point out certain differences between 
this and the corresponding result obtained in Chapter III. The reader 
may recall that the present problem differs from the case treated in 
Chapter III only in that the present problem considers the additional 
requirement that the total power output should match a specified level 


p, at time t,. 


Comparing equation (44) and equation (38) of Chapter III, which 
constitute the optimality conditions for both problems, it can be noticed 
that the additional requirement manifests itself in the form of two 
changes introduced in equation (38). The first is the addition of a new 
term, aK, (t) Ki); to the kernel K(t; a) of equation (III, 38). The 
second effect is the addition of rs Ean TIT,F y to the forcing function 


fee in equation (III, 39). 


The advantages of the present optimization approach over more 
conventional techniques are self evident. Techniques based on variational 
principles and modal expansion methods would yield necessary conditions 
for optimality in the form of an infinite system of ordinary differential 
equations with mixed boundary conditions. In contrast, the necessary 
and sufficient condition (44) constitutes a finite system of integral 
equations, with as many equations as there are control devices in the 


reactor core. 


Equation (44) is amenable to application of successive approximation 
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techniques and, as it was the case in Chapter III, a convergence analysis 


of the contraction mapping algorithm could be carried out with relative 


ease. 


3.5 An Example 


Consider the slab reactor model of the problem example given in 


section 4.3 of Chapter III. 


In the neighborhood of the equilibrium distribution 9 (r)5 
Zee, 
90) = bu 5 sin(*) 
(45) 
the state deviation w(r, t), at any time t > to» is determined by 
v(r, t) = G(r, t3 r', to) Z2,(r') + FG, ts t) UCT) (46) 


where 
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Zy (x) is the state deviation at time to and the eigenvalues 


ga are defined by 


2 
a eae dhe 2 nDV 
2 dail 2 ~ a] a0 ¢ (49) 
Tt D b 


Associated with Z 9 (r) there is a total power output Po» given by, 
Po = &¢ re f Zy (x) dr (50) 


If in addition to the requirement of the problem in section 4.3 


in which the control U(t) minimizes 


ic 
vf 
su) = f fw G, t) drdt + ky fo UM(r) Ulr)dt (51) 
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iE 
0 
the control function is also required to satisfy the power constraint 


b 
te 0 = eZ, a Wir, t,)dr (52) 


+ 
Te and 


then it follows from the previous discussion that the operators T,> l 


x } : 
qT) constitute the only new concepts required to fully characterize the 


necessary and sufficient condition for optimality. The particularization 


of these operators to the present case is given below. 


In view of (52) and (11) it follows that 
1 Ap a 


SEpuat 
T.U = e.25 f f G(r, t.3 re Tt B(r)U (1) drdr 
0 to 


b 
= e,f. Pt F(r, ty3 TIT IAT As (53) 


> iS 
So ‘ + 
feylavisgte ay brs g? Sed 
1 
(€s) 
vd asviga ..9 Juqseo Ts 
(zy — 
Fiat “~ f ; y Sn Tkupsy atts o2 polsihbe: i 
“ = a ‘4 J R 
slntaln (nu foutben 
« 
q 
[ ons ee ’ 6 x o 7 
(it) rier) o ! bib 3 yah 9 L- L* 4 
ee 
am 
an Yawod sid Ee} 3a : % mes prod. avte ek notaswe § 
¢ 4 - i a ’ ” 
(Sz). woes .2w \ 3737" 
A 7 
bal 24 : 
Rey . 
; % 
bras Ty LT #7OUSTITO any ioia notesvosth sieatves mw; ens ao) 
| ee 
od ssisetnerads viv}. oo barlupes. 2dgenaoo oan hao 
« 1 , " ‘ ok a ‘ 


seainaaae ~H#iioitise eae .¥iilamt > oes | sul baee beens 


wolod sevie ef wae snag 


2g sunlion. 2 


Furthermore, making use of (48) one finds that 
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It may be recalled from section 4.3 that B(r) denotes the M- 


dimensional row matrix operator with entries 
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u,(t)dt 
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Given that G(r, t; r' t) is self-adjoint and in view of (36), it follows 
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that the adjoint Ty 
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with entries Ty (t) given by 
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zt (2n-1) 


(56) 


is an M-dimensional,time-dependent column matrix 
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(2n-1)tr, (2m-1)tr, See 
Soiree 21 esi 
2 tan am—1) (1-89) : 
EDO, A p— | oR 
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x *_-] 
=T_[T_T. ] that the pseudo inverse 


Finally, it follows from qT, 1T17, 


“0 Ne ; 5 ; : : oe é 
Ty is a M-dimensional column matrix with entries qT, (t) given by 


1 


ty, 
T, () = = (58) 


It can be shown that the kernel [a °K, (t)K,, (a) - K(t; a)] in 
the optimality condition, for this particular example is an M-dimensional 
matrix with entries Ky .(t3 Pak is a) in the £th-row and jth-column 


given by 
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Computation of the optimal control would require that the kernel 
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in the optimality condition be replaced by an approximate version, 
which for example, could be obtained by truncating the series in the 
defining relations (59) and (62). In addition, the convergence rate 


of the successive approximation algorithm would have to be analyzed. 


Although the present case is more involved, the derivation of 
error estimates between the true and the approximate kernels, and also 
the computation of the convergence rate of the contraction mapping 
algorithm, could be carried out in the same manner as it was done in 


Chapter III. 
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at Some Comments on the Problem of Controlling the State Distribution 
During Load-following 


A natural variation on the previous problem would arise if the 
requirement that the total power output from the reactor core matches 
a specified level at time ty is replaced by the more stringent constraint 
that the total power output should match a given load trajectory at a 
finite number of points, in the time-interval of interest. This 
formulation, clearly, has relevance to the problem of controlling the 
reactor core during load following operations. Although the previous 
and the present problem can be treated on an equal footing by the 
optimization technique of the minimum norm, the latter is, of course, 


more involved computationally. 


It may be recalled from the previous section that the relation 
connecting the control function U and the total power output p(t, ); at 


time ths is given by 
= ; 63 
E(t.) = T(t,s Punt) (63) 
where E(t.) is defined by 


E(t.) = p(t,) - J weet, Be sel Ateneo 1G: 


and T(t, 3 t) defines the transformation from Hy into the real line, 
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It is clear that if the power output is specified at the finite set 


of points {t,t}; BUS? 25) o's Ke? along the time-interval [ty> : then 


1)> 

the function E(t, ) would also be completely specified in {t,}. 
One of the typical features that characterize the optimization 

techniques of functional analysis is that the transformations involved 


in the formulation of the problem are defined on normed function spaces. 


In the present case, the function E(t, )5 bet Ole N. defined 
in (64) can be seen to be an element of a finite-dimensional linear 
space. In fact, it could be considered to be an element of the Hilbert 
space Hy» (an N, -dimensional Euclidean space) with inner product 

Nx 
ia abce hen delineed Chey Bity)on- (66) 


If € denotes the column vector in Hy with entries E(t, ) defined 
by (63) and T represents the column matrix transformation with entries 
T(t, 3 t) defined by (65) then the present control problem could now be 


reformulated in the context of functional analysis. 


Find the control Ue Hy that minimizes 
su) = |[Fu-ylly + [Ul ly (67) 
2 af 
and satisfies 


ee aT Thc (68) 


Again, it follows from Porter's result that the necessary and 
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sufficient conditions for optimality are given by 
* a x 
Ups OF FU + Tan + Fy, (69) 
and 


ane Tana, 
- TF FU+n+TF y (70) 


wy 
i} 


: hae ; 2 
where now n is an element of H, and T is the pseudoinverse of T. 


4 
It should be noticed that in this case the pseudoinverse involves 
A Ak 
the inversion of the N,.~dimensional matrix defined by T T . Where as 


usual, the asterisk denotes the adjoint. 


“+ 
For N,. small, the computation of T poses no problems and the 


conditions (69) and (70) may be combined to yield the operator equation 
Geet Te Lr FU +A Cia) 
in which the function A is defined by 


A=Te - (TT - I) Fy (72) 

However, for large values of Ny the computation of [tr]! 
becomes more difficult and some suitable iterative method would have to 
be used. For this reason it is convenient to rewrite conditions (69) 
and (70) in a form that does not involve, explicitly, the pseudoinverse. 


To this end, the function W in Hy, is defined as follows 


f Da’ 


4. 


P 7 - ‘a 4 
; =) ‘ ian ee, “The 
“T io Serevatnweeq ofa et 7 bus Be Jorg ie Le \ 


iy soyeQgontobpoatc art , obfe a gens best yor a4 bivode | 
¥ 
a cs i 


7S ity ¥ Dood bbrlteh xiviem ieoolene ptidties il oil? Ae. rol ‘ 


intofhs gle asoneb’ oe. 


Ee 


6 eastdorg aq. esepq { 30 nobteneywen aid Elna, 


" LIV pPS: TOI KISGS sHs yt tet ¥ OT? OG ores eal ¥an san) bite (0a) 4 


hae | is 


- @8%) 


oe: bo mantuan age 3 
ad it age Rev i 


iwi 


tiers we : wisn i : 


° 


es. ° au 
i 


4 
of 
ue 


70 


TT W=n (73) 
Substituting (73) in (69) and (70) the desired result is obtained: 


* “k * 
=Sr Og Uptariwit Fry (74) 


G 
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and 


ereeoui tol) TOW + To Fy re (75) 
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Finally, introducing the parameter € in order to rewrite (75) in 


a form that lends itself to application of successive iteration techniques, 
One. nO te * 
Wea se TariraU: — fe .T T .- Ll W+ielé-TF.y] (76) 


in which the parameter e€ may be used to improve the rate of convergence 


of the iterative algorithm employed. 


It is interesting to point out that in the limiting case where Ne 
tends to infinity, the Hilbert space Hy, would no longer remain finite- 


dimensional. The inner product would be defined by 


t 
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est oe = fe eer (cjede (77) 
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and the transformation T would simply become T, as defined in (11). 
The optimality conditions would be given in this case by expressions 


similar to (74) and (76) where instead of T, the transformation T is 
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written. It is not difficult to show that in this limiting case, the 


x 
adjoint T is given by 


* 
TUT; t) cat) = 


t 
— 1! - — 
ky jee * f Br()Q™(r)G Cr, ts rf 1)Q *c'yH (ear E(t)at 
T Vv 


(78) 


This section is ended with one last but important observation. 
The reader should notice that the problem discussed above is based 
on the implicit assumption that the specified load trajectory lies in 
the range of the transformation T. This assumption is of course 


required in order to ensure that the equation 
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ae Optimal Control of the State Distribution with Fixed End State 


Another variation on the problem of controlling the state 
distribution along a given time-interval would arise if in addition to 


the requirement that the control function U should minimize the functional 
5) = |[Fylly + [IUlly as) 
2 nf 


the requirement that the state should reach a specified distribution 
Z, (r) at time ty is also considered. It is clear that this formulation 
would contain the case treated in section 3, in which only the total 
power output is required to match a given power level at the end of the 


time-interval. The present problem, however, would be much more 


restrictive and consequently much more involved computationally. 


The additional restriction imposed on U would take in this case 
the form, 


Z, (x) ri Cree tere Me et ed 


i. 0 = F(r, th3 t) U(t) (80) 
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where all the variables are defined as before. The optimality conditions 


would be given by 
k k * 
U = - F FU + FW sr Sy (81) 


and 
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where Fi denotes the operator F(r, t13 t), mapping Hy into °H, (Wi ie 


an unknown function in H and — is a function defined by 


E(r) = Z, (r) = EK t 3 ig ty) Zo (83) 


x 
The adjoint Fy is given by 


Fi(ts r')w(e') = ko ’R f Bl (rya@(eye" cr, ts tr} t) W(e")dr 
V 


(84) 


and the other variables in (80) and (82) are defined as before. 


Although the conditions (81) and (82) are amenable to the application 
of successive iteration techniques, the degree of complexity posed by 
these conditions is formidable. In order to show how difficult it 
would be to implement an iterative algorithm, it suffices to realize that 
at each iteration step, in addition to the trial control function, at 
least one state distribution defined throughout the reactor core and 
the time interval [ty> t,] would have to be stored in memory. This 
requirement severely limits the application of the present problem 


formulation. 


The degree of computational complexity that characterizes this 
problem could be reduced considerably if the requirement that the state 
of the reactor core should match a specified distribution at time ty 
is replaced by the less restrictive requirement that only a projection 
of the state should match a given distribution in a finite-dimensional 


subspace of H at time ti: Specifically, if the state distribution 


V(r, t,) is represented in terms of a complete basis {py} in H, 
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foe) 


ptr, t,) =a) x, (ty) ©) (85) 
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where {x} is the set of expansion coefficients, then the formulation 


of the less restrictive requirement would take the following mathematical 


form: 
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H 
(86) 
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where {y 3 is the associated-dual set of {ty} satisfying 
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and M, is the dimension of the subspace of H into which the state is 


projected. 


It is clear that in the limit as M, tends to infinity the require- 
ment (86) would become identical to (80). For the particular case, 
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for example, where {py} and {y 3 are chosen to be the natural modes, 
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where A is the operator that appears in the state equation 


st (r,t) = A(r) W(x, t) + B(r) UCt) , (90) 


the coefficients x, ft) of the expansion 


foe) 


v(r, t) = L x(t) v(x) (91) 


would be given by 


A (t-t,) M t A_(t-T) 
n 0 * n * 
x(t) =e <Zo> at fy if e <B,> Dar Ag heen 
j=l to 


(92) 


where BS is the jth-column of B and u, is the ith-entry in the control 
vector U. Also, the linear constraint in (86) would take the form: 
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where the {x, } are the coefficients in the series representation of 
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at): Suboptimal Control 


So far we have not directed our attention to problems in which 
the space of controls is finite dimensional. This situation arises 
naturally when digital computers constitute the main controller in a 
feedback system, or artificially when either time-discretization or 
function-expansion techniques are applied in order to simplify a 


particular optimal control problem. 


In what follows, it is assumed that the control function is given 
by a finite series expansion in terms of a set of known functions of 
time Xt), cont PS AA M, that is 


Ge eo ( th) (94) 
1 | 


where fa, 3 are the expansion coefficients. 


It may be noticed that this concept applies equally well to the 
cases where the control function is either discrete or continuous. 


The particularization to either case depends, obviously, on the nature of 


the functions tw}. 


In view of the modal expansions (91) and (94), it follows that 


the cost function (8), 
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can be rewritten in the form: 


t 


1 
J(U) = - 
(U) ul oy 2 z,(t)) ep, (a(t) - 2, (t))de 


ftp ea be 


+k (95) 


ui a bs ee J “nj Mnjki “ik 


where 26. is the jth-expansion coefficient of the desired distribution 


pir, tyes 9. 18 given by 
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i Sas m 
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where Toi is the ni-th entry in the matrix R. Similarly, substituting 


(94) in (92) one finds that the expansion coefficients x are given by 
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In order to reformulate the optimal control problem in the context 


of functional analysis it is convenient to introduce the following 
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function spaces: He Shall represent the Hilbert space of time dependent 


sequences, with inner product 


1 2 »@ 


re yom fey” ya(e¥no Bx it) ode (99) 
ae 9 n=l j=l J ya Re es 


where the coefficients toe are defined in (96). A. will denote the finite 


dimensional Hilbert space of real matrices with inner product 


ee De } eee) 2 ey leseyere nan (100) 


where the coefficients Byki are defined in (97) and the variables O 
and ik denote the nj-th and ik-th entries in a and $ respectively. 
Finally Hp shall represent the finite dimensional Hilbert space with inner 


product 
© 1c eget ie (101) 


aie The Minimum Norm Formulation 


With these definitions the suboptimal control problem can now be 


formulated in the following form: Find the control element a in He that 
minimizes 
AV) a 
su) = ||Fa-yll, + [loll (102) 
A, H, 


and satisfies the linear constraint 
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where F is a linear transformation from H. into H. with entries defined 


by 
Peat 
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T is a linear transformation from Ha into Hy with entries Fftp3 
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and y is a function defined in H by 
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From Porter's result it follows the unique solution to this problem is 
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% UM — a 
Seite rio [i'n + Fy). (109) 
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The adjoint F is defined by the inner product relation 
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Substituting (104) in the left side of (110) one finds that 
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and introducing a new index p which uniquely corresponds to pairs (j, k) 


it follows that (111) can be rewritten in the form 
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In view of (112), one finds by inspection that the adjoint F is 
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Pe, given in (113). 


: > Vk 
Similarly, it can be shown that the adjoint T is given by 


Vk -l1 -1 
T § = ky B 8€ (115) 
where 8 is a M3 4 M, - matrix with entries 
s) <5 . 
= > 
pn ae be H Ban (ty) (116) 


in which p uniquely corresponds to the pair (j, k). 


Finally, it follows from (114), (113) and (104) that the operator 
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[F F] is a matrix of dimension M, defined by 
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The reader may have noticed that H. is in fact a finite dimensional 


subspace of the Hilbert space Hy» which as it may be recalled is endowed 


with the norm 


JJul[2 = a u(t) R U(t)de (120) 
Hy 0 to ° 


It follows from this observation that if the functions fw, (t)} are 


dense in H 


1 


DP then the optimal control in He is a suboptimal element of 


H, that tends to the optimal control as the number of functions M, in (94) 


tends to infinity. In order to show that this assertion is in fact true, 


it is convenient to point out that minimizing the performance index (8) 


is equivalent to minimizing 
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Given that the functions Wy are dense in H_, for any positive 
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real number o, there is a positive integer M, such that 
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Assuming that F is bounded, it follows from (125) that 
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e=([||Fl|tl]lo. (128) 
If U is the element of He which minimizes (121), then 


3) < 3M) < JW) te, (129) 


and since ¢ can be made arbitrarily small, it follows that J(U) tends 


7 r 

+ (yi ita Bekite | 
: og! Cet aah LB (wt > 
[rr 7 . : ‘ : 

Vio? as e y 

bem 9 te aan wit . en . 

aly estaned {3h 

x ] 


rs 


a 4 


he 
Anis af 
(TF? AREY we larteet . | 
co apeee & i> gae: phe 3 
ov cae 20 
. ; a vs a ‘+ 
Ph 7 ys a%, 
>” ‘ 


84 


to J(U) as M, tends to infinity. 


ee An Example 


Consider again the linear, homogeneous slab reactor model of section 


4.3 of Chapter III, 


2 
Sir, eJucevD ase t) + V[vz,-E.] v(r, t) 
M 
ay a u,(t) d9(r) 6(r-r,) (130) 


with boundary condition 


vw, t) = pdb, t) = 0 (231) 


and the steady state distribution given by 
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Also, for simplicity, Q and R are assumed to be the identity 
matrices. The desired trajectory Z(r, t) and the desired final distribution 


Z, (x) are both chosen to be zero. 
BAT Cc) = Z,(r) = 0. C135) 


If the functions pe (e)s are chosen to be the unit pulses 
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and the index p is given by 
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Also, the matrix F F is defined by 
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where m and k take on integer values in [1, M,] and j and i take on 
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The adjoint T is given by 
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6. Numerical Results 


The slab reactor model considered in the following numerical 
examples is based upon the data given by Wiberg in reference ([11], p. 359). 


The complete set of data used in these examples is presented in Table l. 


In this particular case, the eigenvalues are real. The first 


eigenvalue is close to zero and all the other eigenvalues are negative. 


For the case where no control action is applied, Figure 1 shows 
the trajectory of the flux deviations along the time-interval [0.0, 2.5] 
seconds. At the end of the time-interval all the higher harmonics have 
decayed almost completely and only the fundamental mode remains virtually 


unchanged. 
The flux deviations are normalized by the maximum value ou: 


In all the cases, the control functions considered belong to the 
finite dimensional space spanned by pulse-functions of duration equal 


to 0.25 seconds. 
Case l. 


Figure 2 shows the trajectory of the flux deviations for the case 


where the control action minimizes the performance index. 


16 t) 


b 2 
f wr, t) drat tk, f J u(t) dt (154) 


0 ty i=1 


and in addition satisfies the constraint that the fundamental and the 


first two harmonics should vanish at the end of the time-interval. 


The optimal control is given in Table 2. 
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Migration Length D/z, = 160 on 
Average Neutron Lifetime av = 0.1 sec 


Infinite Homogeneous 


Multiplication Constant vi /E = 1.0256 
Neutrons/Fission v = 2.5 
Diffusion Coefficient D = 0.5070 cm 
Reactor Width b = 250 cm 
Time Interval [tos tj] to = 0.0 sec 

ty = 2.5 sec 
Weighting Parameter ko = 10.000 
Number of Rods M=2 
Number of Linear 
Constraints in Equation (153) M, = 3 
Duration of Pulses w(t) [t,-t)1/M, = 0.25 sec 
Control Rod Location oa 66 cm, ty = 150 cm 
Initial State Distribution Z)(r)= b (0-05, (r)+0.03p, (r)+0.01, (r)] 
Desired State Distribution Z(r) = 0.0 


Table 1 Data For Numerical Examples 
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Figure 1 The Uncontrolled Flux Distribution 
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Figure 2 The Controlled Flux Distribution: 


Case l. 
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Case 2. 


Figure 3 shows the trajectory of the flux deviations for the case 
where the control function is only required to minimize the cost 
functional (154). Although the results of Figures 2 and 3 are almost 
identical, it is clearly evident from Figure 3 that at the end of the 


time-interval the fundamental mode is still present. 
The optimal control for this case is given in Table 2. 


Case 3. 


Figure 4 shows the trajectory of the flux deviation for the case 
where the control is required to satisfy the constraint that the 
fundamental and the two first harmonics should vanish at the end of the 


time-interval and, in addition, should minimize the performance index 


2 
Pac) (155) 


The effect of the optimal control function upon the trajectory 
compares poorly against the first two cases, as should be expected. The 
optimal control function computed in this case constitutes the effortless 
way of satisfying the given constraint on the final state. The control 


function is given in Table 3. 
Case 4. 


Figure 5 shows the flux deviations for the case where the control 


function minimizes the performance index 


t 
yi -12 
JU) = f w, t))dr +k J 2 u, (t)dt. (156) 
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Figure 3 The Controlled Flux Distribution: 
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Figure 4 The Controlled Flux Distribution: 
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Figure 5 The Controlled Flux Distribution: 
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The performance of the optimal control in this case does not compare 
favorably against the first two cases where the deviations of the flux 

are penalized along the trajectory. Although this result should be expected, 
it should be mentioned also that the comparison of the present against 

the first two cases is not entirely valid, since for the same weighting 


parameter k, in the performance indices (154) and (156) the absence of 


0 
the time-integral in the first term of (156) implicitly gives more weight 
to the second term of the performance index, which of course penalizes 

the control effort. Consequently, the norm of the optimal control in this 


case is almost two times smaller than the corresponding norms in the first 


two cases. The optimal control is given in Table 3. 


The computations were performed on the Amdahl 470/6 digital computer 


of the University of Alberta. 
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CHAPTER V 


SYMMETRY REDUCTION OF REACTOR SYSTEMS AND 


ITS APPLICATION TO OPTIMAL CONTROL 
i Introduction 


Although the discussion in the previous chapters concerning the 
application of the minimum norm techniques applies to linear distributed 
parameter models in general, the practical implementation of the derived 
theoretical results to more realistic reactor models, involving for 
example geometrical configurations in 2 or 3 dimensions and including 
the dynamics of fission products such as xenon and iodine, can only be 
made possible if efficient methods for treating both the mathematical models 


and the resulting conditions for optimality are utilised. 


One of the main difficulties that arise in the case of realistic 
3-dimensional reactor models is posed by the large number of control devices 
involved. As the reader may recall from the discussion in the previous 
chapters, the conditions for optimality derived for a variety of control 
problems include at least as many equations as there are control devices 
in the core. It is not uncommon for large reactors to have 14 or more 


of such reactivity control devices [96]. 


This chapter is devoted to the problem of decoupling, or reducing, 
either the mathematical reactor models or the derived conditions for 


optimality into more tractable and independent submodels of lower order. 


The reduction approach presented here takes advantage of the 


geometrical symmetry that may exist in a nuclear reactor configuration. 
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Its foundations are based on the well-developed theory of the represent- 


ations of abstract groups [34, 94, 95]. 
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Die Some Basic Concepts of Group Theory and Symmetry Principles 


The discussion of this section is aimed at establishing the 
operational notation that will be utilized in the chapter. Only the basic 
concepts of group theory that are needed to describe the approach to 
reducing linear distributed parameter systems by means of symmetry 
principles are cited. The detailed discussion of definitions and group’ 


theoretic results is left to the references [34, 94, 95]. 


An abstract group is a set of distinct elements with a binary 


multiplication law such that 


a) Multiplication is defined for every ordered pair of 
elements of the group, is closed, single-valued and associative. 
b) There exists an identity element I in the group. 


c) Every element has an inverse relative to the identity I. 
The number of elements in a group is called the order of the group. 


A representation, I, of a group, g, is a group in which the 
elements are matrices. The matrices in IT are associated with the elements 
of the group g through a many-to-one correspondence which preserves multi- 


plication. 


Two sets of square matrices, say {D,} and {Dj} POreiE ls 2 alae ft 


are said to be equivalent if there exists a square matrix M, such that 
, nis : 
Ds = M D M i=l, fds eeey h (1) 


for some ordering of the matrices {Dit. 


A representation of a group is said to be reducible if it is 
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equivalent to a representation of g composed of matrices in block- 


diagonal form. 


The character of a representation with matrices {D,} is the set 


of quantities {x,} where x, is the trace of the matrix D,. 


Ee D, and D, are elements of a group then the element De D, D, 


is said to be the conjugate of D, with respect to D 


2 i 


A class of a group is a maximal set of mutually conjugate elements. 
The main theorem of the representation of groups states that every 
group has exactly as many inequivalent irreducible representations as 
there are classes in the group. Furthermore, if r 64) = (D5) and 
9 ORE 9 ee Bee 
T = {D, }; j=1, 2, ..., h are any two of these representations, then 


the matrix entries satisfy the following orthogonality relations: 


h 

alee e, DY DSP lm 7 0 AE ate (2) 
PR Cahe hq). 1 

b) oy it [D, Nye eet) ais n#i (3) 


and/or m#k. 


ec) Y vie ppordy 2 2B (4) 
bee j j nm R 
j=l ~-m q 
where p 6 denotes the mn-th entry in the jth matrix of the representation 


mn 
rq) , ea, denotes the inverse of Dae hand te represent the 


order of the group and the dimension of each of the matrices in the 


(q) 


representation [. » respectively. 
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The symmetry of a body is described by giving the set of all those 
transformations which preserve distances and bring the body into coincidence 


with itself. Any such transformation is called a symmetry transformation. 


For a given body a complete list of symmetry transformations satisfy 


the group properties. This is called the symmetry group of the body. 


If a body is finite in extent, the symmetry transformations can 
all be built up from two fundamental types: 
a) Rotation through a definite angle about some axis. 


b) Mirror reflection in a plane. 


In what follows, the ordered collection C = {C, } will represent 
the set of reactor regions, the union of which constitute the reactor 
core. The regions are disjoint and the only common points between any 


two contiguous regions are located at their shared boundaries. 


A symmetry transformation defined on the core, which physically 
represents a positional interchange of reactor regions such that distances 
between points of the core are preserved, can be seen to represent in 
the context of abstract sets a permutation operation defined on the ordered 


set C. 


If R denotes a symmetry transformation on the core then RC = {RC, } 
denotes the image of the set C under R, and it is composed of the same 


collection of regions but now in different order. 


For example, the rectangle of Figure 6, which arbitrarily has been 
decomposed into four equivalent disjoint regions, has its symmetry 


described by the following symmetry transformations: 


a) Oo reflection with respect to the y-axis. 


ha 
b) Ons reflection with respect to the x-axis 
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c) R, rotation of mt radians with respect to the center. 


d) I, the identity. 


Applying the transformation R upon the rectangle, one finds that 


the regions 1 and 3 are interchanged, and so are the regions 4 and 2. 


Therefore, making use of our notation, it follows that, 


RC, = Cy 
RC. = Cy 
RC, = Cc, 
RC, = Cc, 


and the ordered set C = {C Cy} is transformed into 


2 es 


3° Cy» Cy» Cot. 


RC = {C 
The dependency of a function w(r) on the spatial variable r, which 

is defined on the region occupied by the core C, will be explicitly shown 

in the form p~(C). This notation is particularly convenient since it 

facilitates the extension of the concept of a symmetry transformation on 

a body, or on an ordered set C, to the concept of a symmetry transformation 

on functions of the spatial variable r. In this manner, if Ris a 

symmetry transformation belonging to the symmetry group of the body C, 

then by definition, the image R w(C) of a function ~(C) under R is given 


by the function ~ defined over the transformed body, RC. This suggests 


the notation 


R p(C) 4 pc) . (5) 
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If ~(C) is an element of some linear space, then it is clear that R 


defines a linear transformation from the same linear space into itself. 


It is pertinent to mention here that the function defined on the 
reactor core C could either be defined on a continuum, as is the case 
for the flux distribution, or defined at discrete locations, as is the 


case for the control functions. 


For example, the distribution ~(C) defined on the rectangle of 


Figure 6 by 
p(x, y) = cos (~~) sin (PY) : (6) 


where x and y take on values over [-a/2, a/2] and [-b/2, b/2] respectively, 


is transformed by the reflection operation op) into 

o4 W(x, .y) = -cos (—) sin (A) ; (7) 
which in our compact notation becomes 

o, wc) = -¥(C) . (8) 


Similarly, if u, denotes some specific property associated with the 
region C,> in the rectangle of Figure 6, then the vector 


U(C) = coL[u, , Uy» Ugs u,] would be transformed by o, as follows 


3 


o, U(C) = coL[u, » Uz> Uys u,] : (9) 


A collection of functions {f,(c)I, finite or infinite, is said to 
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to generate a representation {D, } of a symmetry group g = {R}, associated 
with a body C, if for every symmetry transformation Ry and every function 


mechs, the term oe, satisfies a relation of the form 


Ree MC eae BOSH eOee co i7C)he (10) 

a Sy - sty k 
After this lengthy introduction of a new notation we are now in a 

position to profit from the results of the main theorem of the represent- 


ation theory of groups. 


Consider a Hilbert space H of functions defined on a region or body 


C which has a symmetry group g = {R, } of order h. 


A representation of g can be generated by applying all the symmetry 
transformations Ry to any function f(C) in H. For example, if f. denotes 
the function £(R,C), then the set {f,1 i=1702. w.0s hideftines e:collection 
of functions that also belong to H. Since for any transformation Ry 
in the group g and any function ; the product Re is also in the set 


{£3 it follows that the product Rife can be represented in the form 


Q 
PeeeeraieeeD.) , f£. |; (11) 
Ce 
where & is the number of unrepeated functions in the set {fj} and ue 


denotes the j-th matrix in the representation of g generated by the 


function f. 


In general, the representation of g generated by an arbitrary function 


is reducible. 


It follows from the main theorem of the representation theory of 
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groups that an arbitrary function f(C) can be decomposed into orthogonal 
components, Paa(C)) which have the property that generate only the nth- 
row of the q-th irreducible representation of the group g. That is, if 
No denotes the number of classes in the group g, and be defines the 
dimension of the matrices in the q-th irreducible representation of g, 


denoted by ae j=1, 2, ..., h. Then the function £(C) can be rewritten 


in the form 


QL N 


q C 
pioeemey. f°) (c), (12) 
n=l q=1 


where the components end k(C) transform under the symmetry operation R, 


according to the form 


Q 
hae eee ey 4 (q)..(q) 
Bee 0) a eu Ft Come (13) 


(q) 


in which I, 
i 
nk 

irreducible representation of g. 


denotes the nk-th entry in the i-th matrix of the q-th 


Furthermore, the components ft (oD can be obtained by applying 


the projection operators 


Ts 
fe ee Ua (TD 1a Row, (14) 
h i i 

=1 nn 


where n and q can take on integer values in [1, at ance (i Nol respectively, 


x 
and the symbol [ ] denotes complex conjugate. 


In summary, the projection operators (pV } decompose the space 
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H into the direct sum of subspaces A : 
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aN, 
fC a he (15) 
n=l q=l 


where the functions in yD 


. generate the nth-row of the qth-irreducible 


representation of g. 
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3 Symmetry Reduction of Nuclear Reactor Models 


The application of the group theoretic results of the previous 
section can drastically simplify the analysis and treatment of problems 
involving linear distributed parameter models which, due to the nature 
of the physical system that these represent, exhibit some degree of 
symmetry. 

Consider for example the general linearized distributed reactor model 


of Chapter III, 


ob (r, t) = A(r) v(x, t) + B(r) U(t) , (16) 


where r is the spatial variable defined on the reactor core C. yp 
denotes the state of the reactor core. For fixed r and t, w(r, t) rep- 
resents an N-dimensional vector. A(r) is a matrix spatial operator 
involving operators of the diffusion type. For fixed time t, U(t) is 

a M-dimensional vector representing the effect of the control rods. B(r) 


is a rectangular matrix spatial operator of the appropriate dimension. 


Also, associated with the equation (16) are boundary conditions of 


the form 

w(r, t) = 0 (17) 
at the boundary of the core, and the initial condition 

v(r, ty) = 25), (18) 


at time te. 
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In terms of our new notation the equation (16) and conditions 


(17) and (18) are given by 


— v(e,t) = Ac) W(c, t) + Ble) U(e, t), (19) 

v(c, t) = 0 at the boundary , (20) 
and 

Wc, to) = Z(c) . (21) 


It should be noticed that although the vector U does not depend 
explicitly on the spatial variable r, clearly, it represents a function 
of the core since every entry in the vector U corresponds to a control 


rod, which in turn is associated with a certain location in the core. 


If the geometrical configuration of the reactor has the symmetry 
group g = {R,} and if the reactor parameters that appear in the mathematical 
model are distributed over the core according to the same symmetry 
described by g, then it follows that the operators A and B are invariant 


to all the symmetry transformations Ro. 


In view of the linearity of the operators A and B it is clear that 


(qd); to 


in this case the application of the projection operators {p. 


equation (19) and conditions (20) and (21) yields 


2 yee, t) = Alo) ¥V(e, t) + Ble) UM (e, t) , (22) 
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VY (ec, t) = 0 at the boundary (23) 
and 
w Mee, t) = 2) (ce), (24) 
n 
where 
MW ua 
od (q) 
ct 1 wy Yn (c, t) (23) 
and 
Cc ac 
U(c, t) = yu’, ¢) (26) 
q=l =p i 


Thus, the projection operators {pi} decompose the space H 
of state distributions and the space E of control vectors into the 


invariant manifolds 6D and nd) 


. L, of A(c) and B(c) respectively. 


It follows from this observation that optimal control problems involving 


a performance index of the type encountered in Chapters III and IV, 


t 
1 
ru) = f tlly-zlly + kollull gla (27) 
t 
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and mathematical models such as equation (16) can be reduced into a 


number of subproblems which consist of minimizing 
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subject to (22), (23) and (24). 


The benefit of the reductive approach is clear. While the original 


problem implies the search for M functions of time, tu, (t)}, in the case 


(q) 


n 


of the subproblems it is the dimension of the subspaces E what dictates 
the number of functions of time that ought to be computed. Depending on 
the degree of symmetry exhibited by the reactor core, these subspaces 
could have a dimension substantially lower than M. 


It should be pointed out that the parameters appearing in the 


(q) 


7. can all be obtained from 


defining relation (14), for the projections P 


the properties of the symmetry group g. 


For most of the symmetry groups of interest in the physical sciences 
these properties have been tabulated and are available in the literature. 


(See for example Reference [14], Appendix IIT). 
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4. An Example: Decoupling of the Optimality Conditions 


The present example illustrates the benefits of using symmetry 
principles in the treatment of the optimality conditions which result 
from application of optimization techniques such as the minimum norm 


formulation. 


As the reader may recall from example 4.3 of Chapter III, the 
necessary and sufficient conditions for optimality that specify the 


solution to the problem of minimizing the performance index 


conor aby a 
Team g(r; tldrdt + ky te) ous (jdt, ) (29) 
to fe) to i=1 


subject to the slab reactor model 
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See ea Pa iy . 
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M 
Ayan CURL) 06 Ce) or cms (30) 
i=1 
bt) =o YE sin) , (31) 


with boundary conditions 


wLO, t)= p(b, t) = 0, (32) 


and initial condition 
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v(x, ty) = Z(r), (33) 


are given by the coupled set of Fredholm's integral equations 


t 
-1 1 es» ee 
u, (t)=A,; (t)-kp f é tae a) uy (a)da, (34) 
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Tene cis 5 OE 
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in which the function KO a) is defined by relations (58) in Chapter 


III. The forcing functions A(t) are given by 
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If the control planes r=r, i=l, 2, ..., M are placed symmetrically 
with respect to the center plane of the reactor; r=b/2, one finds that 
the core would have a symmetry group g composed of two symmetry trans- 
formations: the identity, I, and the mirror reflection with respect to 


the center plane, o. 


That {I, o} satisfies the properties of a group is clear from the 


relations 
II=I1 (39) 
Io=olI=o (40) 
and 
oo=1 , (41) 


This group is identified in mathematical physics as the symmetry 
group c.. (Reference [14]; Chapter III, describes the group classification 
system used by chemists. Although the classification of symmetry groups 
is irrelevant to the present discussion, its use may simplify the process 


of obtaining from compiled tables the properties of a given group). 


The group g has the two unidimensional irreducible representations 


given below 


(q) 


The projection operators PA 


are in this case given by 


= [Ito] (42) 
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and 
p{2) -4 Pci. (43) 


For simplicity, suppose there are only two control rods (planes) 
symmetrically located with respect to the center plane. Then it follows 
that the kernels Ke, in the optimality conditions (34) would satisfy in 
this case the relations 


“a a 


Kyo(ts a) = Ky, (ts a) (44) 


and 


aA 


Ky, (t3 a) = Kyo (ts Ae (45) 


Also, since the effect of applying the transformation o on the reactor 
is felt upon the optimality conditions as the interchange of subscripts 
1 and 2 it follows that the kernels Ky are invariant to o. Therefore, 
the reduction of the optimality conditions into two decoupled equations 


by means of symmetry considerations is possible. 


With only two control functions, u, (t) and u,(t), the space E 


becomes the two-dimensional Euclidean space. Applying the projections 


(q) 


n 


Ee are in this case the unidimensional subspaces spanned by the normalized 
geal yee 
vectors 0 “V9 and P1 Vr respectively. Thus, by expanding 
35 1 -1 


to any base in E one finds that the invariant manifolds ae and 


the control vector U in the form 
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No 


vey = ¥  BSY(e) p64 , (46) 
1 1 
q=l 
where the {a1} are the time-dependent expansion coefficients, and 
substituting (46) in the optimality conditions one finds after taking 


the inner products in E that the conditions for optimality become 


t 
ope Seay a Wied ay tar 
Bog eCe) om B [A, (c)+4, (1) I-ky u [K,, (3 a) 
Ke (0) 186 (ada (47) 
738 Dae 1 
and 
t 
ee (2) 4. Ce) Jka ae (ts a) 
iL /z 1 2 0 ty iB 
- K, (ts 0) 18}7? (a)da (48) 


These equations can be solved independently. Given the simplicity of 
the present example the decoupling of the optimality conditions could 
have been achieved perhaps by other than the sophisticated group theoretic 


techniques. 


However, the reduction of mathematical models in the case of reactors 
with 2 or 3-simensional configurations becomes more difficult and it is in 
this case that the benefits of applying a systematic approach are welcome. 
The next example illustrates the application of the symmetry principles 


to a three-dimensional reactor model. 
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Di. An Example: Reduction of a Cylindrical Reactor Model 


Consider the homogeneous cylindrical reactor of Figure 7 and the 


one-energy group model 


3 
ater, 8, z, t)=VD mee B 2emt eo V[vE_.-2 lv(r, ae 


yi B(r, 0, Zz) Ucey (49) 
where 
14 
B(r, 6, z) U(t)= V oo(r, 8, 2) ) S(z-z,) a, (r, 8) u(t), 
=I 
(50) 
and 
Or, 8, 2) = 4, cas(> z) J,(y,, 2). (51) 


The function a, (r, 8) is equal to 1 in the small region occupied 
by the ith neutron-control device and equal to zero elsewhere. v 
denotes the Laplacian operator in cylindrical coordinates, L denotes 
the length of the core. Yo, represents the first root of the zero-order 
Bessel's function JQ (YR) » and R is the radius of the cylinder. All the 


other variables have the usual meaning. 


The flux is equal to zero at the boundary of the core and the 


initial condition is given by 


pr, 0, Z, ty) = Z)(rs 8, zie (52) 
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Given the symmetrical allocation of the control devices in the 
core, as shown in Figure 7, one may conclude that the symmetry principles 
discussed in the previous sections can be applied to the problem of 


reducing the mathematical model. 


The reactor core has the following group of symmetry operations: 
the identity I, anticlockwise rotation of m1 radians with respect to the 


x, y and z-axis; denoted respectively by Ri» ie and Ro 


Reflection with respect to xy, yz and zx-planes; denoted respectively 
by Se oes and Fx And finally, the inversion operation with respect 
to the center, Ri: 
These symmetry transformations, together with the group multiplication 
table shown in Table 4, constitute the symmetry group Day of mathematical 
physics. Each of the symmetry operations forms a class in itself. The 


eight unidimensional irreducible representations of D are given for 


2h 
reference purposes in Table 5. It follows from this observation that the 


mathematical model of the core can be reduced into eight subsystems by 


the projection operators oes sad Oped ARE PON ars 


In this particular example, the space E of control vectors becomes 
the 14-dimensional Euclidean space. The vector U = cok[u,, ree 14! 
transforms under the symmetry operations as shown in Table 6. For 
simplicity only the subscript indices are shown. The Table shows for 
example that the 6-th and 11-th entries in U are interchanged by the 
operation Ro 


(q) 


By applying the projection operators PL to any complete set of 


(q)) 


vectors in E, one finds that the eight invariant subspaces {E; 


q=1, 2, ..., 8 are spanned by the vectors coed?) shown in Table 7. 
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Table 4. Group Multiplication Table. Sample: R R =F 
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Table 5. Irreducible Representations of D 
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Also, the eigenfunctions of the spatial operator 


4 
VDV" + V[vi--2] (56) 


generate all the irreducible representations of the group Dane Table 8 
shows the classification of the eigenfunctions according to the irreducible 


representation which they generate. The set of eigenfunctions is 


orthogonal in the real Hilbert space H, endowed with the inner product 


TD Te 42 
<p, o> = f f. f xv v@,e, z) o(r,6, z)drdedz. (57) 
220i .0 


By applying the method of eigenfunction expansions one can find that the 
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invariant subspace Hy 
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spanned by the infinite set {py » is given by 
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where QD denotes the infinite set of subscripts associated with the 


eigenfunctions in eye and we denotes the eigenvalue associated with 
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Laplacian Modes w 
ni 


£ Cos ((2n-1) Fae) cos (18) J, Y5,F) 


J T a 
f ,Sin(2n L z) Sin (ié) J, (4,0) 


T z . 
f Cos ((2n-1) 7 z) Sin(ié) J, (V5.7) 


; 1 : 
f,,,Sin(2n z) Cos (i 6) J (V5.8) 


- T ; 
f ,, Sin(2n z) Cos (ié) J, (V5.0) 


T aera? 

f Cos ((2n-1) L» Sin(ié) J, (¥4,7) 
: 1 aes 

f ., Sin(2n L 2 Sin (i6) J, ¥4,4) 


f Cos ((2n-1) a Fe Cos (i8) J, (V5.8) 


Where 
ee bore, 35 Kis Gee a3, 
-2<2<¢ 0< 6 < 2n 
Vat is the kth root of J, (yR) 
f 2/(vmL R rel (7;,2)) if i#0 
i sg 
Y2 | (VnL R Tsay (75,2) if i=0 
(q) 


Table 8 The Laplacian Eigenfunctions Voik aes 
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GE es T 
Lee (2n-1) L° LOr-d#l_ 3940.55 (60) 
or 
(q) 2nt 
Be Me fOr, =2. a5 (61) 
a denotes the expansion coefficient corresponding to the base- 
vector pv) in the expansion of the control U, 
N 
8 q 
Deemer) 6.2) (x) pc (62) 
m 1m 
q=l m=1 
and Ne is the dimension of the invariant subspace Ewe 


While the dimension of the control space E in the original model 


(1) (4) 
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CHAPTER VI 
A MODAL EXPANSION APPROACH TO LOAD-FOLLOWING 


i Introduction 


Given the large difference in time constants between neutron 
kinetics and reactor poisoning dynamics, the analysis of problems such 
as xenon spatial stability or load-following can be simplified 
considerably by making use of the valid assumption that the neutron 
kinetics is an instantaneous process, compared to the slow dynamics of 


the core poisoning effects. 


With this assumption a very general reactor core model would 
consist of the multigroup neutron diffusion equations at steady state, 


coupled with the xenon and iodine dynamic equations. 


A new method is presented in this chapter, for obtaining an 
approximate solution to such a model. The objective is to obtain an 
explicit functional relation between the neutron flux and the neutron 
absorption cross sections of the zone controllers. The motive is to 
reduce the number of variables involved in the optimal control studies 


of load following operations. 


The proposed approach is described in terms of the one-energy 
group model with xenon and iodine dynamics. Several examples are 
presented in which use is made of data that corresponds to a typical, 


pressurised tube, heavy water-moderated, large nuclear reactor. 


134 


WIOS-TADI OF 
i 
— 
. h tg 
; yajad efaetevo: 9nid BE Bo eee ad? t 
“us 
ams} »%y Shevtans odd ,cotmengbhs antHoetaq SoIeRet | ben's 


in 
i! 


tle 9 en gniwoilot-heot 5 ysiirdase inna > 2 
tg- 4003 ne Pique : bifsv 984 We Bap gutaan ef ids nab aoe 


yeryrimiinitens 6 et be 


A) 

~ 

mM 
] 


; *j/1 ; ou 
aghetees gutnoekog 
i we 
I, . fayberm eto 5" x Vrev & Rolo giomes aids a2. 
He (6h oon set gic so hslum nts Yo 
Sy i. an 
meAare ex rbor brie CHAK abs nate 
ape syaiqur: stil ek begeotanty St bodtoer i 


. Wes 
se: ! sido oa ,fsbea 6 dows Ad potswlos S7amkze 

’ es, ; : 1 ‘sf ai 
: oft eawied soltaie® benofsandt Jieke 


rorsued @/17 pre iil: _ os ; i oH 
: ol hee 

i ad *i ° i + ie) ee 
6) at evitom aif  .erellowtaos sak Say 29 snottose @s07: aotte 


: ioe 


anttura tae fnos tpakaes ais al boytlownmt nalivanam be redea 76 
al ; ‘ 


aor ; enabrerage 8 


¥ 
‘ . | : 
aig asiTAxSa [a3eyse 25 lite ao A b itieatl nth rato 


tuokgvi & 01, aie agrarian teas eonk Ye aha 
“i 

; r eat iy’ 

bn ctasebee soe 


Sérh-oh. sel Be eased mt f Sab saith et doorways ae qd 


SHiosss, “— egret 
= uA 7 i-w y 


Ue he 


De Background 


The variety of reactor models used in connection with the 
control of distributed nuclear reactors are usually classified as either 
nodal [1,2,3], finite-difference [4,5,6] or modal expansion [7,8,9] 
depending on which method is used to obtain an approximate description 


of the reactor core processes. 


Often a method classified as one of the above is further categorized 
depending on the details and specific features of the approach. A 
modal approach could be categorized for example by the type of modal 
expansion used. A description and comparison of the basic approaches 


to distributive reactor modeling is given in reference [10]. 


The nodal and modal expansion methods are, by far, the most widely 
employed approaches to reactor control problems in which use is made of 
distributed reactor models [11-20]. The method of finite-differences 


is used extensively and almost exclusively in reactor physics calculations. 


The main advantages of the modal expansion approach over the 
finite-difference and nodal methods can be listed as follows: 

a) The modal expansion method is capable of representing the 
distributed nature of the reactor core in the same detail as would be 
obtained with the finite-difference approximation. The number of equations 
required is, however, much less for the modal expansion approach. 

b) The number of equations required to represent the distributed 
reactor is comparable for the modal expansion and nodal methods. The 
spatial representation of the reactor core is however more detailed for 


the modal expansion approach. 
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In general, the main disadvantages of the modal approach are the 
following: 

a) For many modal expansion methods it is not possible to 
determine or even estimate error bounds between the true solution to 
the problem and the solution obtained by means of the modal expansion. 

b) For some modal expansion methods the spatial modes are 


difficult to compute [10,15]. 


Several methods which make use of explicitly known, analytical 
functions as expansion modes have been proposed in the past to circumvent 
these disadvantages. The method of solution functions of R. Bobone 
[21-23] is probably the best example. The method was used successfully 
in criticality calculations.and also in computing the steady state 
power distribution in a reactor core. His approach assumes the homogeneity 
of the diffusion coefficient and reactor buckling within each of several 
core regions. The method consists of expanding the neutron flux within 
each region in terms of solutions to the diffusion equation and evaluating 
the expansion coefficients so as to minimize the mean squared error of 
the neutron flux and current at the region interfaces. This method, 
unfortunately, cannot be applied to the analysis of load following 
operations, since the xenon distribution within the reactor regions makes 
the neutron absorption cross sections not only spatially but also temporally 
dependent. Iwazumi and Koga [18] proposed a method in which the spatially 
dependent parameters of the diffusion equation are replaced with constant 
terms that are chosen to keep the fundamental mode steady. The spatially 
dependent terms were then treated as the initial distribution of the 
control variable. The basic idea was to use the eigenfunctions of the 
Helmholtz operator as expansion modes. This approach was used by 


Iwazumi and Koga in connection with the problem of changing the flux 
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distribution from a given initial state to a desired final state in a 
short period of time while minimizing a functional that penalizes the 
terminal flux distribution error and the deviations of the control 
variable from the initial control distribution. The one-energy group 
model with one delayed neutron precursor in a slab reactor was considered 
in their work. Recently, Iwazumi's method was used by Asatani et. al. 
[19] in connection with the same control problem and the same reactor 


model. 


Iwazumi's method has not been extended to treat the load following 
or xenon oscillation problems. Whether or not this extension would 


lead to justifiable simplifications is an open question. 
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Ss The Approach 


The method presented in this chapter can be classified as a 
modal expansion approximation in terms of analytical functions and it is 


characterized by two important features: 


1) The use of analytical functions as expansion modes makes 
it possible to bypass the problem of computing the expansion modes by 
numerical means. 

2) The formalism of functional analysis, which is used in 
the development of the approximate model, makes it possible to estimate 


bounds on the approximating errors. 


Three basic assumptions are considered: 
1) The geometrical configuration of the reactor core is such 


that a complete set of wavefunctions in L,[V] is explicitly known and 


2 
has the property that the functions vanish at the extrapolated boundary 
of the core. L,[V] is the space of squared integrable functions over the 
core volume V. 

2) The neutron diffusion coefficients for the different neutron 
energies are assumed to be homogeneous throughout the core. This 
assumption is acceptable for large reactor cores that in addition to 
being lightly rodded either do not make use of reflectors [36] or have 
reflectors made of the same material employed as neutron moderator [36]. 
The multigroup bucklings may be spatially dependent. 

3) The reactor model is linearized in the neighborhood of a 


given power distribution which could either be at steady state or also 


be a slow transient during a large power level change. 


The description of the proposed approach is given in terms of the 
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one-group model with xenon and iodine dynamics presented in Table 9. 
For the sake of clarity, the procedure has been divided into five steps: 

1) The neutron diffusion equation (1) is transformed into a 
nonhomogeneous Fredholm's integral operator equation of the second kind 
[24-29]. 

2) The kernel of the resulting integral operator is approximated 
by a degenerate kernel. An error bound on the L,[V] norm is obtained, 
of the error between the true solution to the Fredholm's equation and the 
solution obtained by means of the degenerate kernel approximation. 

3) The integral equation is solved for the neutron flux by 
applying the method of degenerate kernels [24-29]. An explicit relation 
for the neutron flux is obtained in terms of the xenon concentration 
and the control variables. 

4) The explicit expression for the neutron flux is substituted 
in equations (5) and (6) and an approximate solution to these equations 
is obtained by expanding the xenon and iodine distributions in terms of 
a finite number of wavefunctions. A bound is found on the Lo{V] norm 
of the error between the true solution to equations (5) and (6) and the 
solution obtained by means of the modal approximation. The bound is 
obtained by showing that the error satisfies a nonhomogeneous Volterra's 
equation of the second kind in the product space Lo[V] and then using 
the Neumann series representation for the error. 

>) Using the explicit expression for the neutron flux and the 
modal expansion solution to equations (5) and (6), the neutron flux is 


expressed in terms of the control variables only. 
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Table 9. One-Energy Neutronic Model with Xenon and Iodine Dynamics 


The model describes the dynamic behavior of the neutron flux 
deviation Vy and the deviation of xenon and iodine concentrations 


from their corresponding equilibrium distributions bo» Xp and To: 


The parameters o, Lo» A. and Vy denote, respectively: the neutron 


I 


absorption microscopic cross section of xenon-135, the decay 


constants of xenon and iodine, and the iodine yield per fission. 


u denotes the control-function deviations from the equilibrium 


controisu:. 
0 
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Vy, (e, t)+ b(E)¥,( t)= E(x, t) 
b(r)= [vz (r)- Mn)S u(x)- ox (x)]/D 
Ee, t= Fo ©, t+ a (r), x, t) 
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Initial Conditions: 


ae (cy ty) sa ¥,0 (r) 
Boundary condition at the core extrapolated boundary (38V): 


Vy Es t) = 0 
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ae Solving the Diffusion Equation 


The one-group neutron diffusion equation (1) can be seen to be a 
nonhomogeneous Helmholtz equation in the Hilbert space LIV] with a 
space-dependent parameter b(r), which by assumption is a piecewise 
continuous function of the spatial variable r. An approximate solution 
to equation (1) could be obtained by expanding the neutron flux in terms 
of a finite number of functions. In fact this procedure is used 
extensively in the literature. The simplicity of this approach however, 
has to be weighted against the disadvantage of lacking error estimates. 
An alternate approach consists of transforming the equation (1) into an 
equivalent integral equation that is amenable to both error analysis 


and computation. 


The first step toward transforming (1) is to express b(r) as the 


sum of two terms, 
b(x)= bot b(r) (10) 


where bo is spatially independent. There are of course an infinite 


number of different ways of doing this. The considerations that dictate 
the best choice will be discussed later. For the moment it is assumed 


that b in addition to being a real number, is also different from all 


th 


the eigenvalues of the following boundary value problem in Lj [Vv]: 
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¥, = antes reav (12) 


where V is the region occupied by the reactor core, dV is the extra- 
FAs 
polated boundary and V’ is the Laplacian operator. By assumption the 


set of functions ty} is explicitly known. 


Using (10) in equation (1) one finds that 
Vv Ce, t)+ by V4 {x C)=er et b (x) V4 E> t) (13) 


and using the Green's function G(r, r') corresponding to the left side 


of equation (13) the following expression is obtained: 


vy t= f 6G, x") £@')av'- f 6G, x") eS bets av" 
V V 
(14) 


where the integrals are over the core volume. 


By invoking the first assumption of section 3 it is not difficult 


to show that the Green's function can be obtained in series form; 


teh Aled Gy 
Cur eria=" ) (pen) (15) 
j Oo Pe 


where {a5} and ty) form the eigensolution to equation (11). 


Equation (14) is a Fredholm's integral equation of the second kind 
in the function space L,[V]. Also, since the integral operator in (14) 
ta compact, it is clear that there exists a unique solution in L,[V] 
provided that the number (-1) is not in the eigenvalue spectrum of the 


integral operator [38]. 
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4.1 The Method of Degenerate Kernels 


Of all the methods available for solving the nonhomogeneous 
Fredholm's integral equation, the method of degenerate kernels [24-29] is 


particularly attractive to use when the Green's function (15) can be 


adequately approximated by a finite series. 


Substituting the Green's function (15) in the integral equation 


(14) the resulting expression is 


v. (x) : 
= Aa 1 ' i ' 1 1 
¥,@ ¢) ta) Pdap EEE t)- bee" )b Ce", t))dv"] 


(16) 

and can be rewritten in the form; 
p. (x) : 
v,(z, t)= : oo a) m) [e, (t)- C(t] (17) 


where the coefficients e, (t) and ae are defined as follows; 
c.(t)= f v,(r) d(x) v,(r, t) av (18) 
j (oil) ae 5? 
e,(t)= f w.(r) f(r, t)dV (19) 
j } i 


and satisfy the algebraic system of equations; 
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J 
where : 
Biryob tr) p,(r)dV¥ 
Riphctpesd seein cn Se (21) 


ij (bot A,) 


a Pu bith y wary 


he tv, 
eo 


ashy Sahn 


eyoensgagnainoa an? gitt vies eid 
-¢ [Of-30) elebrsd atactarapst ic oe 


ai nae CEL) qokdon? # baer 23 


Fares | asta 
fe y ns 
osideune’ isvesgnk edy mt (Ex) cork sam Nga 


sawed. ton. 62 hi, Ta me (hy? tiethe kan ta as 
6 i | ' 


rt ve la aye wwe i 


a i 


(a) haf “4 | “G7, 


Rohl taupe ee 


A natural approach to solving equation (16) consists of approx- 
imating the Green's function by a finite series and then solving the 


algebraic equation (20). It is clear that the approximation can be 
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made as accurate as desired by just including more terms in the series. 


In addition to the computational simplicity of this approach the method 


also has the property of being amenable to error analysis. 


The objective in this section is to determine an error bound 


between the true and the approximate solutions to the equation (16). 


AV) 
It is convenient to introduce the following notation; G(r, r') 


denotes the approximate Green's function and 2 denotes the finite set of 


subscripts {j} corresponding to the terms that appear in the series 


nu 
representation of G(r, r').. That is; 


p. (r) ps fe ) 


A (22) 
(byt x) 


av) 


avy “a a 
Aisa, 1190 1s a function in L,{V], then L, L, L and L denote the 


following transformations in L,{V]5 


Lo 4 f G(r, r')o(e')av' (23) 
V 
to A f Gr, x')o(x')av' (24) 
V 
to A ubo = f G(r, r')b(e')o(z' av! (25) 
V 
and 
Pe ene w 4 
Lo 4 tbo = f G(r, r')b(x')$(x")dv' (26) 
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version are written in operator form; 


¥,fe)= et (t)< 1 ¥, (t) (27) 


and 


av) 
vy (t)= L £(t)- L W(t) (28) 


% 
where Vy Y% and f are shown as functions of t to emphasize their 


temporal dependency. 


From equations (27) and (28) the error e(t)= Ves s(t) is 


given by the expression; 
av ry st 
Set =e L-L jer (t)— [L eer L V4 ft) ] (29) 


which can be rewritten as follows 


AV) 


e(t)= (L-L) £(t)- L e(t)- CL), (t) (30) 
Hence the L,[V] norm of e(t) is bounded above by 
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and it follows that if L is a contraction operator, then the error 


estimate is obtained in the form 
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Given that the set ty} is complete in the Hilbert space L(V], 


every function ¢ in L,{V] has a fourier series representation 


= d Wass ¢ > (33) 


where are o> denotes the inner product in L,[V]; 
(x) $(r)dv (34) 


and ¢ is the complex conjugate of ¢. Also, using the inner product 


notation it can be seen that 


v.sv., £(t)> 
Aci, 00 Bes Beat tee 


(L-L) £(t)= J (35) 
jE Sere 


and given the orthonormality of the functions Oe 


(L-L) £(t)= , (E(t)- 2 ¥,<d,, F(e)>)> (36) 
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= (L-L) (£(t)- J vy<vy, £(t)>) (37) 
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Similarly using the relations L=Lb and L=Lb, it follows that 
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Oe ey Abe oes iis 
(L-L)V4(t)= (L-L) (bY, (E)- 2 Vy<by» bY, (C)?) (38) 


ieN 
Substituting (37) and (38) in the inequality (32), the error estimate 


is finally obtained: 
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| |L-] | 
[le(t)|| < —— get) (39) 
[1-| |u| [1 


where 


g(t)= |l£(t)- <y., £(t)> tees) at oe 
| | £ (t) a ViSY, ||+] [by, ce) Duke bv, (t)>| | 
(40) 


The norms of the operators L, L, if and (L-L) are given by; 


Plas ax, (41) 
j [bot «| 
||| |= max —+— (42) 
jen [bot | 
(ical ks 3) eee (43) 
j#2 |by-r, | 
and 
[|| |= max y,7/? (44) 
j 


where ty;} is the eigenvalue spectrum of L*L and L*¥ is the adjoint of L. 


“a 


Often it is difficult to evaluate the norm of L and only an upper bound 


estimate can be obtained. A conservative estimate is given by; 


[|z|| < max —+—— max |b(r)| (45) 


j [bot | xr 


The derivation of these norms is presented in Appendix II. 
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Inequality (39) provides a useful criterion to determine how many 
and which terms should be used to approximate the Green's function (15) 
by a series expansion. It is also clear from the same inequality (39) 
that bo» in the decomposition of b as shown in (10), should be chosen 


so as to minimize the norm of L. 
4.2 An Example 


The data for this example has been obtained by applying the modified- 
one-group approximation to the two-energy group data presented in the 
Appendix I, The modified-one-group data and the constants corresponding 
to the xenon and iodine dynamics are given in Table 10 . The core 
configuration and the layout of zone controllers considered here is the 
same as in the example in Chapter V. This is shown in Figure 7. 

The core is cylindrical and exhibits a degree of symmetry that corresponds 


to the point group D of mathematical physics [34], which has eight 


2h 
unidimensional, irreducible representations. 

The Laplacian modes for the cylindrical case {bak are ordered 
with the help of three subscripts. These modes are presented in Table 8 
where these have been classified according to the irreducible represent- 
ation of D,.. The set of eigenvalues 12 ak? corresponding to the 


2h 


Laplacian modes {Yak are given in Table 11. 


The parameter b, was chosen to be the constant value required to 
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Migration area 
Diffusion coefficient 


Infinite multiplication 
factor 


Xenon decay constant 
Iodine decay constant 
Iodine yield constant 


Absorption macroscopic cross 
section of xenon 


Neutron flux distribution 
at steady state 


Core configuration 
Length of core 
Radius of core 


PIFSte root 1oL Jy (YR) 


Table 10 
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D= 1.2349 cm 
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neighborhood of small radius and located in the center of the core is 
negligible then it can be seen that for most flux shapes of practical 


interest the estimate (49) is indeed valid. 
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x The Modal Expansion 


From the discussion of the previous section it is clear that the 


neutron flux LF can be expressed in terms of the function f in the 


operator form; 


gp aad 8 (60) 


where F is the operator (or an approximation to) [I + host In the 
particular case of the example, the equation (60) corresponds to the 


equation (56). 


Substituting equation (3) in (60) one finds 
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and eliminating the flux from equations (5) and (6) the dynamic equations 


for xenon and iodine are obtained as follows, 
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where 
De coL[y, v1 (65) 
z= coL[z,, z_] (66) 
ee geet (67) 
ae = x oFu (68) 
and 
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I Iniep °° 
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At this stage it is convenient to pause and make the following 
remark; note that the linear analysis of xenon spatial oscillation could 
be carried out by computing the eigenvalues of the operator A. In fact 
this procedure would be equivalent to studying the xenon stability by 
means of Kaplan's natural modes [7]. Several difficult problems, however, 
of theoretical and computational nature make this approach impractical. A 
simple computation would show that the operator A is not compact, non- 
self adjoint and furthermore, not even normal. The theoretical implication 
of this, is that nothing can be said about the existence and completeness 
of the eigenvalues of the operatorA. Even if it is assumed that the 
eigenvalue spectrum of A is discrete, the numerical computation of the 
eigenfunctions and their corresponding eigenvalues is difficult. Successful 


computations have been reported only in the case of a slab reactor SEL 
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Because of this a number of methods have been proposed to obtain approx- 
imations to the natural modes. Of these, probably the best known are 
the Lamda and Mu methods developed by Stacey [10,15,30,31]. Roughly 
speaking, these methods consist of Laplace transforming equation (64) 
and manipulating the resulting equations in a way that allows the use of 


conventional, statie diffusion codes. 


The objective in this section is not to perform a stability 
analysis but to find an approximate solution to equation (64) in terms 
of the control variable u. It is obvious however that the two problems 


are interrelated. 


The method proposed consists of expanding the iodine and xenon 
concentrations in terms of a finite number of the Laplacian modes 
ty} and obtaining a weak solution to equations (63) and (64) in the form 
of an integral equation. Let I denote the finite set of subscripts 
corresponding to the modes in the finite expansion and let v, and v. denote 


the approximate solutions, 
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and 


d 2 

—W (t)= A.W. (t)- Fs oO ee $2 x F 

dt Xs I I, cle v, ( ee ee D *o eke es 
+ <Wey 2, (t)> (73) 


with the intial conditions; 
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In order to evaluate an estimate of the approximating error it is 
convenient to combine (76) and (77) in an operator equation similar to 


(64). To this end the following definitions are introduced; 
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the form; 
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a Eveaa yt 2 (84) 
where 
AV) 
r A 
a if if 
A= (85) 
av) 
eer Ay 
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the equation (86) could be written in the form 
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where 
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h = (A-A) b+ (2-2) (90) 


Then using the familiar convolution integral, 
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‘Equation (91) has the familiar form of the nonhomogeneous Volterra's 
equation of the second kind in the product space L,[V] x L,[V]. tis 
well known that a unique solution to (91) exists provided that the 
operator o,(tst)A is bounded in L(V] x L,[V] for all finite values of 
(t-t). Furthermore, the solution can be expressed in the Neumann series 


(contraction mapping series); 
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K, (t,t) = o Ct, tA (94) 

and 
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and by induction; 
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which tends to zero as n tends to infinity. Using (100) it follows 


from (93) that 
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Even though this error estimate is very conservative, together 
with (90) and (92) it can be of help in determining how good the approx- 
imation v is. In particular, (90) is a good criterion for checking 
the adequacy of the approximation. Note also that the sharpest estimate 


obtainable by this approach would be determined by choosing A, se as to 


minimize the bound M in (96). 


5.1 The Functional Relation 
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(105) 
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where Arp Ary» Ayy and Ayy are square matrices of dimension equal to 


the number of subscripts in QI and where 
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recalling from (61) that 
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Hor 5 eave a (61) 
and since py = ) vs W. » the neutron flux is obtained in terms of the 
LeQl <f 


control variable u; 


eee 5 Fe.) We (eet mule) (111) 
LeQr i 


where for simplicity only the temporal dependency is explicitly shown. 


Substituting (110) in (111), the desired functional relation is 


obtained; 
eee kote a | 6+ Fu (112) 


where A is a function in L(V] defined by 


A(t) = TEE io Ciese mea re ran Tle Ca se. 
or Bea eas 4 dy i 20 
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and F. is the transformation from the space of control functions into 
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5.2 An Example 


For the particular case of the example in section 3.2 the matrices 


asd. ,Q) , (a) (q) 
Art : Aly ‘ t and Ayy corresponding to the q-th subsystem are diagonal 


and two-dimensional. The dia 1 i i (q) (q). 
gonal entries in Aly and Ayy 4 


ae (q) (q) 


at) . Raitt . Ge f : oF vaike (115) 
and 
(q) cy D a (q) 
a = A + + — q 
XK (nik) Ratk » 6 es has D me Eos) Vaik if aes 
are given in Table 17, in which the eigenvalues ate and pag of 
(nik) (nik) 
the augmented matrix 
(q) (q) 
Ar Ary 
CLLZD 
(q) (q) 
I 


are also shown. 


It can be seen from the eigenvalues of Table 17 that although all 
the subsystems exhibit oscillatory behavior, all but one of the subsystems 
(4) , (6) 
are stable. The mode Yi01 is unstable. Of all the stable modes, Vii 


exhibits the slowest transient response with a time constant of 47.13 


hours. 


The neutron flux Y% is obtained in the form 


y yD (118) 
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where 


Ne By Sd, ed) ~ (ayy ata (119) 


Expanding the control variable ia in terms of the symmetry-adapted 


linear combinations {eV} k=1, oe of Chapter V3; 


C 


q 
=1 


(q) 


and using the spectral representation of F as given in Table 15 


the terms a p 6D 5d) and Faeaey are given as follows; 
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also 
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where the coefficients are given in Tables 17, 18 and 19. 
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where the coefficients are given inTables 17, 18 and 19. 
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Table 18 
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Number of Rods 
Rod Diameter 
Rod Length 


Location of Rod Centers 


Rod 1 150 cm = 0.6747 Rad 
Rod 2 150 cm = 7/2 Rad 


Rod 3 150 cm 1/2 Rad 


Rods 4 to 14 are located symmetrically as shown in Figure 7. 
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CHAPTER VII 


CONCLUDING REMARKS 


1 Summary 


The optimization technique of the minimum norm formulation in 
Hilbert spaces has been applied to a variety of control problems in 
which a linear distributed parameter model is utilised to describe the 


dynamic behavior of a nuclear reactor core. 


All the control problems considered here are variations on the 
same theme; the control of the reactor's state distribution in the 


neighborhood of a given state. 


To the best of the author's knowledge, only the solution to the 
problem of minimizing a performance index that penalizes both the 
deviations of the reactor's state and the control effort has been 


published in the literature by other authors. 


While the more conventional techniques approach this problem 
through modal expansion methods and invoke variational principles which 
yield necessary conditions for optimality in the form of an infinite 
system of ordinary differential equations with mixed boundary conditions, 
the minimum norm formulation employed here yields necessary and sufficient 
conditions for optimality in the form of a finite set of coupled integral 


equations which can be solved by means of iterative techniques. 


It has been shown in Chapter IV how the additional constraints, 
imposed on both the total power generated in the core and the terminal 
state, can be treated by the optimization technique of the minimum 


norm. Also, it has been demonstrated how these additional constraints 
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can be related so that the resulting problems are suboptimal. 


The application of the minimum norm approach to all the problems 
considered in the thesis was described in terms of a general linear 
distributed parameter model. Several examples were given in which the 
particular case of a homogeneous slab reactor with a finite number of 


controls and with the one-energy neutronic model is considered. 


It was demonstrated in Chapter V how by merely studying the geometrical 
symmetry of a nuclear reactor core it is possible to reduce the 
mathematical model of the core into more tractable submodels of lower 
order. This simplification of the problem will facilitate the practical 
implementation of the minimum norm techniques to more realistic reactor 
models which may involve, for example, geometrical configurations 


in two or three dimensions. 


A new method was developed in Chapter VI for obtaining an approximate 
solution to a reactor core model that is suitable for describing the 


dynamic behavior of the reactor core during load following operations. 


The method, which can be classified as a modal expansion approach 
using the eigenfunctions of the Laplacian operator, yields the solution 
to the model in anoperator form that is amenable to application of the 


minimum norm technique. 


The proposed method has a scope of application that covers those 
cases where the reactor core is of simple geometric configuration, so 
that the eigenfunctions of the Laplacian operator are explicitly known, 


and where the diffusion coefficients are nearly homogeneous. 


The method was described in terms of the modified one-energy 
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neutronic model with xenon and iodine dynamics. 


Both, the symmetry reduction procedure and the modal expansion 
approach were combined in an example that involves a cylindrical reactor 
and utilizes the numerical data typical of a large nuclear power reactor 


of the thermal type. 
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z Suggestions for Further Work 


The application of the minimum norm technique to situations in which 
the saturation of the control functions is also taken into account 
constitutes a natural extension to the present work. While in theory 
this type of problem can be formulated and treated in the context of 
functional analysis, the implementation of the optimizing algorithms is 


difficult. More work is needed in this area. 


An ideal approach to the optimal control of linear distributed 
parameter systems should combine the computational advantages of the 
methods of functional analysis with the engineering benefits of the feedback 
control concept, which characterizes the Kalman-Pontryagin theory of 


optimal control for lumped linear systems. 


It would be worthwhile to investigate, in this spirit, the 
feasibility of synthesizing practical feedback control schemes for 


distributed parameter systems via functional analysis. 


Also needed is a sensitivity analysis of the optimal control 
algorithms reported here with respect to model inaccuracies and system's 


noise. 


Further work is needed to evaluate trade-offs in distributed model 
sophistication versus worth and practical implementability of the optimizing 


algorithms. 


A different variety of problems, that is also of interest to the 
electric power industry, to which the application of the minimum norm 
is worth investigating arises in the area of fuel consumption optimization 


in continuously fueled nuclear reactors [98, 99, 100, 101]. 
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The basic problem in this area consists of finding the optimum 
power distribution, and the corresponding neutron flux, throughout the 


reactor core so that the fuel burn up is maximized. 


Although variational calculus has been applied to some particular 


cases in this area, it may be possible that the minimum norm technique 


leads to more accurate and efficient algorithms. 


The short-range economic operation optimization of electric power 
systems with nuclear, hydro and thermal generating plants is another 


promising area for application of the minimum norm technique. 
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APPENDIX I 


Nuclear Data for a Typical 1200 MW(th) Natural Uranium, Heavy Water- 


Moderated, Pressurized-Tube Nuclear Reactor. 
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Infinite Multiplication Factor 


Resonance Scape Probability 
Neutron Slowing Down Length 


Thermal Diffusion Length 


Fast Neutron Group Diffusion Coeff. 


Slow Neutron Group Diffusion Coeff. 


Fast Neutron Speed 

Slow Neutron Speed 

Neutrons Produced Per Fission 
Core Radius 


Core Length 


Table 21 Two-Energy Neutronic Data 


* Clean Reactor Data Provided to the Author in Private Communication 
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1.404027 
0.89801 


134.3 


235.4 
1.3643 
1.2349 
10/ 
3x10° 
2.640 
400 


600 


cm 
cm/sec 


cm/sec 


cm 


cm 


by Dr. G.M. Frescura of the Ontario Hydro's Department of Nuclear 


Studies and Safety. 
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Fission Macroscopic Cross Section de “wae 
Absorption Macroscopic Cross Section ea = Dth/L? 
Absorption Macroscopic Cross Section ae = igi, 
Scattering Cross Section L =D eon iia 
1/2 f s 


Table 22 Nuclear Cross Sections for the Two-Energy 


Neutronic Model 


Delayed Neutron Relative Abundance Precursor Decay 
Precursor B,/8 Constant A 

4 -1 i 

7 (sec ) 

il 0.0405 0.000763 

2 0.2129 0.0305 

3 0.1355 0.1180 

4 0.3882 Oo: 

5 0.1268 15226 

6 0.0461 MN he 


cee 2 B. = 0.004982 
i 


Table 23 Delayed NeutronPrecursor Data 
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APPENDIX II 


Derivation of the Norms of L L 


Us (Lely and L. 


For convenience the Dyadic notation is used. That is, the symbol 


o> denotes a function ¢ in L(V]. A functional wp in L,[V] is represented 


in the form <j. 


The spectral representation of L is given by 


=f <u, 1 ae 
L Pa ai v5 ( ) 


where 


rE = (Lis 2) 


Also the Fourier representation of $> is of the form; 


= - ae > CLES 3) 
o> P Maas ae p 
= Lee 
t vic ee ( ) 
J 
where 
a ; (Fie) 
we < a9 ) ? 


Given the orthonormality of the functions V5) it follows that 


=I). “ake (II.6) 
Lo> = 2 eee S, j 


j , 
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and 
1/2 
[zo] | = £2 (e2c%)] cErens 
j Sted 
< Max le, | lol . (II.8) 
j 


Note also that the norm of L is defined by 
[Iz] | = sup Il || - (II.9) 
||| [=2 


Consider a function $ with Fourier coefficients C.> where; 
Max |€, 
x (6, | 


ieee ee. se 
Bs i Ue (II.10) 
C, = 4 


max 


and ‘Uso is the subscript corresponding to the term Max 5, I> eat he re Re pe 
i 


le. | = Max |e,|. (II.11) 
JInax * | J 


It is clear that Bev Hise also 


| |Lo|| = Max Esl. (37¢125 
j 


In view of (II.9) and (11.12) it follows that 


[It] | = Max 6, | (II.13) 
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similarly, 


~ 
||u]| = Max |e.] . Ga kaiy 
jen 


A A ~ 


Consider now the operator L = Lb. The norm of Ld> is given by, 


<i, Lg>t/? 


| [L¢| | 


2 (ITV15) 


<p, L*LO>" 
where the adjoint L* is obtained in the form 


L* = b L (II.16) 


a “a a &M 


L is compact, Lb and L*L are also compact. Since L*¥L is self adjoint 
it follows from the spectral theorem for compact, self adjoint operators 


that the spectrum of L*L consists of a finite or countably infinite set of 
eigenvalues and that the sequence of eigenfunctions of L*L is orthogonal 


aA A 


and complete in L(V], provided that zero is not in the spectrum of L*L. 
Therefore, if {U3 and ty,} denote respectively the sequence of eigen- 


functions and the spectrum of L*L, then each function ¢ in L,{v] has the 


Fourier representation 


o> = 2 Og aus p> Gm ign FF A) 
j 


“a “aA 


Also, the spectral representation of L*L is given by; 


L¥Ine SU .> y.<U, (II.18) 
1 Ug fee 
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therefore, 


1/2 


| |Lo|| = <o, L*Lo> 


Deer <Uss b>clans 
1k ee 


Bere ll 


[A 


Max Ys 
j 
and it follows that 


Wiha tee ve 1/2 
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CLIGTS 5 


(II.20) 
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